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INTRODUCTION 


The  purpose  of  the  computer  code  described  in  this  volume 
is  to  find  the  solutions  of  the  eigenvalue  problem  given  bv 
the  differential  equation 


^Up(3,f)  f CO*-  / i _ 7 _ 

l cVz(f)  ^ )]'ur(-z'P)sC> 


( 1. la) 


and  the  boundary  conditions 


(O,  p ) - O 
- Of 


(1.1b) 


(1.1c) 
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The  variable  is  considered  a parameter.  The  solution 

consists  of  sets  of  eigenfunctions 

.(■*,  ? ) jo-1,2,...  N If) 

and  their  associated  eigenvalues  kp(^)  for  each  desired  value  . 

Inputs  are  the  sound  freouency  and  a function  c(  2 p ) 

^ 27?  ' \ 

corresponding  in  this  case  to  the  measured  or  simulated  sound 
velocity  distribution  in  some  part  of  the  ocean. 

In  the  program  it  is  assumed  that  c(  p ) is  given  in  the 
form  of  a two-dimensional  grid  of  values : 

ft.  ) , --  C-(21/( C) 


and  that  one  can  interpolate  linearlv  between  grid  noints.  The 
set  of  values 


FP  , 


-<‘  * i 
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will  be  referred  to  as  the  sound  velocity  profile  at  range 

■ -r) 

point  py  . In  practice  series  of  sound  profiles  are  usually 

. — 9 

taken  along  a straight  line  course  so  the  vectors  p.  can  be 
replaced  by  the  scalar  range  variable  py  , 

The  problem  as  outlined  above  arises  as  oart  of  the 
attempt  of  solving  the  wave  equation 

)V.  Vp  fr>]  * - 'S  (r~ro) 

describing  the  pressure  filed  of  a unit  point  source  (located  at 
/<o  ) in  an  inhomogeneous  ocean  medium  of  density  p ft"  J and  a 

propagation  constant  k(r)  = oj/^c(t)  where  c(r)  is  assumed  to  denend 
arbitrarily  on  the  depth  z and  arbitrarily,  but  gradually,  on 
the  horizontal  range  coordinate  p=(X,^). 

Assuming  the  density  r)  to  be  constant  throughout 
horizontal  layers  and  defining  a velocity  potential  <£>  (■*-*)  such  that 

f(r)=  f Cr  ) 2 4 (;>2) 

one  obtains  for  each  layer  a Helmholtz  equation 

VL<f>  CT  ) + kZ(P)  4>  (r  ) J /P~Pa). 

\ JL  * <j  / 

For  a range  independent  sound  velocity  distribution,  i.e.  for 

-9  * 

c(r>  = c(z)  the  solution  of  (1.3)  can  be  written  as  a normal 
mode  sum  of  the  form 

<f?  Cr) 

r (i.4) 

The  "range  functions"  <\|y,  ( ip1)  then  satisfy 

[Vf%  kf  3 fp-po)  ^Pf?o)y  Q 5) 

which  is  easily  solved,  while  the  "depth  functions"  Up  fid  are 
obtained  by  solving  (1.1)  with  p=ro*.st.  For  arbitrary  but 
gradual  range  dependence  one  can  attempt  a solution  in  the  form 
of  an  almost  senarated 


(1.3) 


(1.4) 


(1.5) 
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normal  mode  sum  which  takes  the  form 


<P  (O  = E Ye  <? 


(1.5) 


The  equation  satisfied  by  the  range  functions  becomes  now 
much  more  complicated  involving  terms  which  couple  all  the 
modes  together,  whereas  the  depth  functions  can  still  be 
obtained  from  (.ti)  except  that  they  are  now  no  longer  range 
independent,  but  gradually  range-dependent  "local"  depth  functions. 

The  problem  of  arbitrary  z-dependence , and  arbitrary  but 
gradual  range  dependence  of  c(r)  is  the  case  for  which  the  pro- 
grams described  in  this  publication  have  been  developed. 

Vol.  1 is  restricted  to  the  discussion  of  the  depth  functions 
(^j?)  and  a method  for  calculating  them.  A discussion  of  the 
general  range  equation  and  a method  for  obtaining  the  nodal 
range  functions  (<p)  together  with  a description  of  the 
associated  computer  code  may  be  found  in  Vol,  2. 

The  method  for  solving  (1.1)  which  was  used  in  the  computer 
code  described  here  is  based  upon  the  assumption  of  oiecewise 
linear  sound  velocity  profiles.  This  assumption  does  not  re- 
strict the  range  of  applicability  of  the  code  in  any  way.  On 
the  contrary,  it  enables  one  to  obtain  analytic  solutions  for 
the  depth  functions  and  thus  makes  it  possible  to  obtain  numer- 
ical solutions  much  more  efficiently  than  by  the  mpthod  of  numer- 
ical integration  which  is  necessary  for  arbitrary  profiles.  In 
this  way,  normal  mode  theorv  becomes  accessible  to  a much  larger 
class  of  problems  than  otherwise  possible.  This  is  of  particular 


-i 

J 
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importance  in  view  of  the  fact  that  the  normal  node  treatment 
is  the  most  exact  method  of  investigation  sound  Dropagation 
problems,  which  still  works  in  cases  where  other  methods 
(ray  tracing,  PE-method,  etc.)  fail.  On  the  other  hand,  by 
inserting  additional  segments  any  velocity  profile  can,  at 
least  in  principle,  be  approximated  to  any  desired  degree  of 
accuracy  by  the  present  method.  This  would,  of  course,  increase 
the  computational  effort  again  to  unrealistic  magnitudes. 

However,  in  practice,  a very  large  number  of  segments  is  neither 
necessary  nor  is  it  in  general  possible,  since  for  any  given  range 
point  there  is  usually  only  a relatively  small  number  of  data 
points  available  to  define  the  local  sound  velocity  distribution. 


L 


II.  DEPTH  FUNCTIONS  FOR  A PIECEWISE  LINEAR  SOUND  VELOCITY  PROFILE 


Since  the  boundary  conditions  are  given  at  z = 0 and  z = cO , 
the  depth  functions  u (V,  » * ) have  to  be  found  within  these  limits. 
The  vertical  sound  velocity  distribution  at  some  particular  range 
point  p which  determines  the  general  solution  of  the  deoth 
function  equation  at  this  point  is,  as  mentioned  earlier,  usually 


approximated  by  a piecewise  linear  function: 

-C  ( "3.  ) = 'CCz:)  t-  S.e  (2-Zt) 


2 ^ 


with 


S-c  =[;c(z;r±)-'C(izl)]/(2^± U-i'ij..  Ktl) 


(2.1a! 


(2.1b: 


The  range  0 4 z 4 oc  is  thus  subdivided  into  horizontal  layers  with 
boundaries  at  z = z^z^for  the  ith  layer.  The  method  to  proceed 
with  is  therefore  to  find  the  solutions  of  the  depth  function 
equation  within  each  layer  and  match  the  solutions  at  the 


layer  boundaries.  The  conditions  to  be  satisfied  are  that 
fj  (2 ) -Up  ( i , P ) 9Axr  p ) /<2 2 


(2.2) 


are  continuous  across  each  boundary.  These  follow  from  the  require- 


ment that  the  pressure 


and  the  velocity 


p ■*= 

Tj  - s=  — \7  (p 


be  continuous  evervwhere. 

The  water  extends  from  z = 0 to  z = z^.  The  water  density 
is  assumed  to  be  constant  everywhere  i.e. 


< (V 


C - £ - ^ -p. 


(2.3) 


The  ocean  floor  is  assumed  to  consist  of  isovelocity  layers 
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with  the  density  constant  within  each  layer.  At  the  present 
time  only  one  layer  is  considered,  hence  the  ocean  floor  is 
characterized  by  2 numbers 

pfe)  = yOQ  'C(%)=  ^ > Z-fe  ' 

However,  if  the  need  arises  these  assumptions  can  easily  be 
generalized  to  several  layers  with  linearly  varying  profiles. 

Within  the  water  c(z,p  ) at  some  particular  p varies 

min  > max(  r ) » where  always  cmax(  /°  ) ^ Cg(  p ) . 

A solution  of  the  eigenvalue  problem  (1.1)  is  obtained  only 
if  the  local  vertical  wave  number  Kp(Z/p  ), 

Kpte,f5)=  {[oj/'cf  M (p  )}  (2.4:; 

is  real  for  at  least  part  of  the  range  0 < z<  zQ( p ) . This  limits 
the  range  of  possible  eigenvalues  to 

Cp  ) < (f)  . (2.5«! 

Discrete  eigenvalues  are  obtained  for 


4 (P  ) > UJ  /c&  (f) 


For 


(2.5t 


4 fp)  < Go/'Cfc  fp), 

any  value  of  kp  will  provide  a solution  to  eigenvalue  problem  (1.1) 
This  report  is  concerned  only  with  the  calculation  of  d'iscrete 
eigeftsolutions . With  the  assumptions  made  above  the  solution 
within  the  ocean  floor,  i.e.  for  the  range 

fp)  4 2 ^ “ 


which  satisfies  the  boundary  conditons  at  ©o  is  simply 


with 


XX  p = ^ fp;  -fr  (p)*} 

Xr  <p'>  - \ ly  - u*/-c|  . 


(2.6a) 

(2.6b) 


- 
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This  leaves  only  the  solutions  for  the  range 

O < * ^ (p  >/ 

i.e.  /Up  (2;p  ),  the  solutions  within  the  water,  to  be  determined. 

They  are  obtained  by  solving  Cl. la)  with  c(z)  as  given  by 
(2.1)  or  more  precisely  with  piecewise  linear  distribution  of 


W'Cf*)]1  : 

ko  te")  = s4  (^~^<)  2^,',! 

(2.7a) 


with 


(2.7b) 


Assuming  a linear  dependence  of  on  z is  compatible  to  first 
order  with  the  assumption  of  a linear  z-dependence  of  c if  the 
variations  in  c(z)  are  small  compared  to  the  magnitude  of  c(z): 

kJVz)  = SioAwi2-  - to'/fsctei)  + S^(Z-^i)}  ^ 
which  with 


A<C-  = I'C  (Zi  + ±)  - * (=£.^  )l  (2.8) 

can  be  written  as 

ko*  (a)  = [u3/-C(Xt-)}a  - 2 

c 2 

Hence,  if  c(z)  is  a linear  function  with  slope  s-  then  k (z)  is 

i o 

also  approximately  a linear  function  with  slope  s.  = -2u>1sf  /{/CCz-i)  V 
and  vice  versa,  if  (2.8)  holds.  Condition  (2,8)  is  satisfied  for 
sound  velocity  profiles  since  typically 

J\  "C  <£.  SO  -m /sec  c = /SoO  -r^/s&c  . 

; 

2 

Proceeding  with  the  assumption  that  k^(z)  can  be  linearized 
according  to  (2.8)  the  depth  function  equation  for  the  ith  layer 


at  some  particular  p becomes 

dxjuf (*, p )/=*»•*■  ■*- { [vA  (*4 / f )'i 2 + s.  (p)('£-'Zi)-^p  ip) j;  u.r p ) = o 

or 

dxMf  (p ; + s«  r/?;  ■*  5 ^ = ° 

with 

dpi  <p)  = fi  (?)*;  - V r^;  . <2.31 

The  substitution 


(*,p  )>  -ftw^rA i-v.-  fp^ (?>*■  *=  wpfe-fccfi] 

(2.10) 

transforms  this  equation  into  the  Airy  equation 
d^/j^7-  = g sUp  . 

Therefore,  the  solution  for -the  pth  mode  in  the  ith  layer  in  the 
water  becomes 

ur?f.,p>  * ■*,<  Ai  (tfi ) - ir;  B,  ( 4>i)  • 

There  are  2 constants  to  be  determined  for  each  layer.  By  watch- 
ing the  solutions  and  their  derivatives  at  the  layer  boundaries 
the  coefficients  of  one  layer  can  be  expressed  in  terms  of  quanti- 
ties from  the  preceding  layer. 


Abbreviating  ) 


A;  = Ai  (?j,i  , 

one  can  write  the  matching  conditions  at  the  boundary  between 
the  ith  and  the  (i-l)st  layer  as 

= Oyi  A*  * **»<  f*<) 


s - - «v[c^  A/  + '1 5 u/" 

Solving  these  two  equations  for  a^  and  b^  in  terms  of  ^(z^) 


(2.12) 
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L 


and  , (z.)  yields 
i-l  i 

£*pi  = n\ulp^±  f* !-)■&/+  ^ u^-.t  f 

"lopi  = 7T  f-^U^  (^,-)Al  ~ \ . 


(2.13) 


The  first  derivatives  of  the  depth  functions  are  obtained  from 
(2.11)  as 


- o/i~iAj.  „* ' ) - -ofi  At  '(Xr ) - V Bt ' r^), (2> 


14) 


keeping  in  mind  that  the  derivatives  of  the  Airy  functions  are 
with  respect  to  % . When  matching  derivatives  at  the  layer 
boundaries  one  has  o-f  course  to  use  derivatives  with  respect  to  z, 
as  was  done  in  (2.12).  However,  in  order  to  propagate  the 
solutions  for  each  mode  from  one  layer  to  the  next,  one  can,  as 
suggested  by  eqs.  (2.13)  and  (2.14),  use  just  as  well  ^-derivatives. 
This  approach  was  taken  in  the  program  described  here.  Defining 

(*■"%  = - <*'“•  V <*’■)  » 


the  coefficients  (2.13)  then  become: 

= ^ < )B/  - T Up/.T. 

V = 1rC -“fi-i  H<)A;'4Tu;,,t  no^  A<] 

T = <*.--./<*;  ■ 

An  alternate  form  for  propagating  the  solutions , also  used  in 
the  program,  is 


(2.13) 


w iih. 


u 


^ u,f)  = '-c±x  ) -t  ^ 


WyK  (%l?\  — Up.c-1  fci  ) 'Cl?-  AXpi-*. 


(2.15) 


wH 


e ire. 
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= -t5  - 

B-c 

h?r)Ai'] 

= rtT{  Bi  (^;  )A^  - 

A* 

C?r; ) ‘B,  } 

^•2i 

%,  )A;'  3 

= fT  - 

- A< 

: '<?*  J i . 

Using  eqs . (2.11^  (2.13)  and  £2 

.14),  the  full 

in  the  water  can  be  determined  in  terms  of  just  two  unknown  con- 
stants, which  can  be  e,g.  a^  and  b , i.e.  the  coefficients 
layer  nearest  the  surface.  One  constraint  for  these 
coefficients  comes  from  the  boundary  conditions  at  the  surfaces 

Up  (0/  p ) = o (2.16) 

Hence  from  (2.11): 


4 4 + ~ O - 

If  one  chooses  for  apl  the  arbitrary  but  convenient  value 

— *^1.  f 

then  b n is  determined: 

Pi 


(2.17) 


(2.18) 


S’4"  ’ (2.19) 

and  one  obtains  for  the  slope  of  the  solution  just  below  the 
surface  from  (2.14) 

U " ' (0,  ? ) = - «i.  (XV-  a = <4  (p)/TT  . 

' ’ * ' (2.20) 

This  choice  of  coefficients  ensures  that  all  solutions  start 

out  with  a positive  derivative  near  the  surface: 

p ) >0  , 

(2.21) 

which  makes  the  X derivative 

(°/?h  ~ U = -TC'1-  . 

r 0 (2.22) 
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Using  (2.11)  and  (2.14)  with  (2.18)  and  (2.19)  one  can  determine 

the  solution  and  the  derivative  at  the  lower  boundary  of  the 

first  layer.  In  terms  of  these  quantities  the  coefficients 

for  the  second  layer  can  be  found  using  (2.13).  In  this  way  the 
w — * 

solution  u ( , p ) can  be  propagated  to  z = z^,  the  ocean 

floor.  There  the  functions  u (^/  P ) and  u ( ^ . p ) have  to 

pi  P ' I 

be  matched.  This  leads  to  the  conditions 


Pw  * pR  up 


(2.23) 


_ YJ  * 

04p  (i/p  ) 


d-^tEL  f). 


(2.24) 


As  indicated  in  (2 -7a) , u contains  one  undetermined 

P 

coefficient.  It  can  be  found  by  using  (2,23).  The  choice  (2.18)  was 
made  because  it  provides  a convenient  starting  point  for  numer- 
ical calculations,  as  will  become  apoarent  below.  But  it  was 
arbitrary  in  the  sense  that  solutions  calculated  with  this  choice 
do  in  general  not  satisfy  the  normalization  condition  for  normal 
mode  solutions,  in  this  case:  \ 


Jo  ^ ^ ~ ^ 


(2.25) 


9 
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To  remedy  this  situation  the  integral  over  the  square  of  the 
unnormalized  solution  is  evaluated  to  find  the  necessary 


normalization  factor 


a„aL  *■  KB  1/" 


(2.26) 


where 


fo  lujf.  {%/P^I  5T  Opt  fit*  (fyt  )J  dfa 

>,•<*«,  r> 


(2.27a) 


= _5 
/W 


p-  e • «• 


27b) 


The  first  integral  can  also  be  carried  out  analytically  using 
-Ta. 

a;  or)  "d?  (Ad?))*]"5*  (J-. 


(2.28a) 


Similarly 

/ A*  (%)3;(z)dg= 


[ ? A«  ) B;  ) - A;'  (? ) Be ' (?  )j 


UltU 

f , , , (2.28b) 

j U fa  (*«a-  fa  WJ  . 

Si. 


(2.28c) 
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These  expressions  follow  from 


f£jt)(^<Z  )*Z  . >^tz)  - CJ  (x)c'J(z) 


by  differentiating  with  respect  to  £ \ 


C,C,  --  C,Cl  (c,‘c^c,c*‘)-C,  'V-C,V."- 


Upon  rewriting  this  equation  as 

cj  lc/'-zc,)  - c/ 

it  is  apparent  that  it  is  satisfied  for  and  being  any 
solution  of  the  Airy  equation. 

-4 

The  functions  u^(  p ) as  discussed  so  far  still  do  not 
satisfy  eq.  (2.24).  The  coeffic  nts  of  the  two  partial  solutions 
u”(  2,  p ) and  u^(  2.,  p ) were  chosen  such  that  the  boundary  con- 
ditions at  2:  = 0 and  t =oo  were  satisfied  and  that  Up(  2? y p ) 

is  continuous  everywhere.  The  derivative  u?(  ) is. at  this 

P 

point  continuous  everywhere  inside  the  water  and  inside  the 

ocean  floor.  However,  so  far  no  provisions  have  been  made  to 

* 

make  it  also  continuous  at  z = z . In  order  to  achieve  this 

B / 

w B 

u and  u have  to  be  adjusted  such  that  (2.24)  is  satisfied.  This 
P P 

can  be  done  by  choosing  the  proper  values  for  kp(  p ) which 
w B 

enter  u and  u through  (2.9)  and  (2.6b).  In  this  sense  eq. 

P P -ih*. 

(2.24)  provides  the  characteristic  equation  for /Eigenvalue  problem 
at  hand.  A more  convenient  form  of  the  characteristic  equation 


14 


is  found  by  dividing  (2.24)  by  (2.23),  which  corresponds  to 
matching  log^arithmic  derivatives  at  the  ocean  floors 

2 Up  fevp)  _ ?»/  w _ v 

9*  “ (0B  Yp  ( P ) MP  C2.29) 

The  characteristic  equation  defining  the  eigenvalues  ( "p  ) 
can,  of  course,  be  established  at  the  boundary  between  any 
two  layers.  It  turns  out  that  sometimes  it  is  advantageous  to 
propagate  the  solution  from  infinity  up  to  the  boundary  between 
two  water  layers.  The  characteristic  equation  is  then 

Up*  pt  (i/i  )U*T  , 

(2.30) 


assuming  the  match  occurs  at  the  boundary  between  the  (M-l)st 

and  the  Mth  layer.  Here  u . is  the  solution  started  at  z = 0 , 

P * 

Up  the  solution  started  at  oo  . 

Both  (2.23)  and  (2.30)  are  used  in  the  paper  and  will  be 
referred  to  later  on.  x 


- > 
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III.  ARGUMENTS  OF  THE  AIRY  FUNCTIONS  AND  GENERAL 
PROPERTIES  OF  THE  EIGENVALUE  SPECTRUM 


Since  the  Airy  functions  play  a central  role  in  the 
problem  discussed  in  this  report  it  is  useful  to  point  out 
some  general  properties  of  these  functions  and  to  try  to  get 
an  idea  of  how  the  various  parameters  which  enter  the  problem 
affect  the  regions  in  which  the  Airy  functions  are  needed. 

It  turns  out  that  the  arguments  of  the  Airy  functions  within 
a particular  layer  can  be  expressed  to  a good  approximation 
in  a relatively  simple  way  in  terms  of  the  frequency,  the 
slope  of  the  velocity  profile  and  a quantity  with  the  dimension 
of  a sound  velocity  which  indicates  how  far  the  eigenvalue  is 
above  or  below  the  wave  number  corresponding  to  the  local 
sound  speed.  This  formula  can  not  only  be  used  to  determine 
the  range  over  which  the  arguments  of  the  Airy  functions  vary 
for  a particular  profile  but  also  to  estimate  the  total  number 
of  modes  possible  and  the  average  separation  between  successive 
eigenvalues  as  a function  of  mode  number. 


The  Airy  functions  are  characterized  by  the  observation 
that  they  behave  like  an  oscillating  function  for  negative 
arguments  and  like  an  exponentially  rising  or  falling  function 
for  positive  arguments.  The  characteristics  are  clearly 
brought  out  by  the  asymDtotic  expressions  for  large  arguments 


( IS  l >>  1 ): 
At  (£  ) - 

Bi  (O- 


c . IS  1 ) 

' 1 ? +,“  t 

g*  r . #sr  / l 

1 1 + 
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Ai(-ir) * — \ H - . i +...)  sin  (y*  ~ ) 

1 S &C*  C V 2*2I6X1*  ' 17 

- ( if  y-")0’5^^  (3.1) 

+ ( i§bk~ " ')*'"■  (v*?)  5 ‘ 

The  argument  of  the  Airy  functions  is  given  by  (2.10)  with 

d .(£)  and  s.(€>  defined  by  (2.9)  and  (2.7b).  It  follows 
P1  » 1 J 

that  (ignoring  the  ? -dependence  and  writing  c(z^)  as  etc.): 


Kt\.  (»)  = - fc)'Vl  [ VI 


(3.2) 


denoting  the  k2  - value  as  a function  of  z along  the  profile. 
As  indicated  by  (2.5a)  and  (2.5b)  the  range  of  discrete  eigen- 
values is  restricted  to 


J^L  - L <k^:k  = — 

'Cu  /ShC, 


This  is  indicated  in  Fig.  1 where  schematically  a typical  profile 


is  shown. 
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The  quantity  which  essentially  determines  the  argument  of  the 
Airy  functions  as  a function  of  z for  a particular  mode  p is  the 


local  vertical  wave  number  K (z)  defined  bv 

P y 


-a. 

P 


(3.3) 


i. 


Defining  new  variables 

( a ) = - + ‘-V'C*,,- 

with  the  ranges 


0 ^ Tc/fa)  $ 


to 


60' 


60  ‘ 


60 ' 


and 


0 (jp*!) 


CO  ' 


to' 


co 


Co- 


a.  2ujx (Aj+Aj,) 


/C  a. 


one  can  see  that 


k/  fa  } = Kf  - - Kj-  fa) 

varies  between 

- 2CO--A*  /-c.,! 

for  p = 1 (i.e.  for  the  eigenvalue  with  the  lowest  mode  number) 
and  at  z = zB  (the  z-value  with  the  highest  sound  velocity),  and 


1 

1 
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for  low  wave  numbers,  negative  arguments  for  high  wave  numbers. 
The  slope  of  the  profile  in  a particular  layer  does  not  affect 
the  sign  of  the  argument.  However,  the  absolute  value  of  the 
argument  may  grow  very  large  as  the  slope  approaches  zero. 
Layers  with  small  or  vanishing  vertical  sound  speed  gradients 
are  frequently  encountered  near  the  surface  and  of  course  near 
sound  speed  minima.  Some  care  is  necessary  to  maintain  the 


accuracy  of  the  solutions  in  these  cases  as  will  be  shown  below. 


In  order  to  get  some  numerical  values  the  following 
parameters  were  chosen: 

*f  « ISO 

= A*.  /&j,  - 0.  i-OO  »»/s«c  loOvvu 

=■  = 5~'m./sec  f>a<  100 

S3  - “ IS**  /sec  per  lOO-m. 

* 4-5*00  •m./se-C. 


The  three  values  chosen  for  s^  correspond  to  low,  average  and 
high  sound  velocity  gradients  as  usually  encountered  at  small 
or  intermediate  depths.  The  value  chosen  for  s°  is  typical  for 
greater  depths.  With  these  numbers  are  obtained  for  the  Airy 
function  arguments  : 

<*)  - -ii.4  Ala),  - o.Sii  a u) , - 0.  dl*  A (2) 

— ZPz  (*)- -o.au  a(*) 

-1-87  A (*)  . 
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For  f = 25  Hz  all  coefficients  are  smaller  by  a factor  of 
2/3 

10  = 4.64.  For  the  lowest  eigenvalue: 

A (o)  s?  — (A*.  + ) 

A ) S - Az 

A ( *0  SO  (3.7a) 

A - Aj  . 

For  the  highest  eigenvalue : 

Alol  »A)  + Ak-Ai-Aj 
A fij)*  A,  + At  - Aj. 

. , (3. 

A A,+Ai 

A (2b)s*  At  . 

The  parameters  encountered  in  a deep  ocean  case  may  be  e.g. 

± s AoO  tw  y “ ^o0  'Tyv  / “ SooO'***. 


(3.7b) 


A^O-H-m/sec,  A^Oon/sec,  Aj-'Ws«,  Ab=3o<w/s«c 


(3.8) 


This  leads  to  the  following  Airy  function  arguments  for  the 
lowest  mode: 


/>*  346  at  the  surface  (z  = 0) 

340  at  the  end  of  the  first  layer  (z  = z^) 
e*  25  at  the  beginning  of  the  second  layer  (z 


0 at  the  end  of  the  second  layer  (z  = z_) 

w 
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0 at  the  beginning  of  the  third  layer  (z  = z ) 

O 

~ 84  at  the  ocean  floor  (z  = Zg) 

The  corresponding  numbers  for  the  highest  mode  are 

-3570,  -3575,  -263,  -288,  -644,  -560  for  the  quoted  value 

of 

-165,  -170,  -4,  -38,  -84,  0 forA^  = 0. 

It  is  apparent  that  a combination  of  high  frequency,  small 
slopes  in  the  profile,  very  low  and  very  high  mode  numbers  can 
result  in  very  large  positive  or  negative  arguments  for  the 
Airy  functions.  From  the  asymptotic  expressions  for  the  Airy 
functions  one  can  see  that  large  positive  arguments  are 
particularly  serious  since  they  can  lead  to  exponents  of  the 
order  of  several  thousand.  For  the  lowest  mode  e.g.,the  Airy 
functions  are  of  the  order  of  unity  near  z = , but  in  the 

example  given, of  the  order 

e**  a io 

near  z = zfc , i.e.  the  eigenfunctions  for  the  lowest  mode  grow 
over  more  than  200  orders  of  magnitude  between  the  ocean  floor 
and  the  location  of  the  sound  speed  minimum.  For  very  small 
slopes  which  may  occur  near  the  surface  or  near  a sound  speed 
minimum  the  exponents  may  be  even  higher.  As  the  small  slopes 

* 

k . — 


do  not  usually  persist  over  large  distances  the  growths  of 
the  Airy  functions  between  the  beginning  and  the  end  of  a 
layer  are  relatively  small,  but  the  coefficients  obtained 
when  matching  the  solutions  to  the  other  layers  become 
small  to  the  same  order  as  the  Airy  function  grows  large^ 
and  near-cancellation  of  very  large  numbers  occurs.  Special 
precautions  have  to  be  taken  to  avoid  overflows  and  under- 
flows in  the  computer  and  to  ensure  that  the  accuracy  of  the 
results  does  not  deteriorate  due  to  the  near  cancellation  of 
large  numbers.  This  problem  will  be  addressed  in  the  next 
section. 


Eq.  (3.6)  can  be  used  to  estimate  the  number  of  discrete 
modes  available  for  a particular  profile  and  frequency. 

The  highest  mode  is  an  oscillatory  function  all  the 
way  from  the  surface  to  the  ocean  floor  since  the  arguments 
of  the  Airy  functions  are  negative  throughout.  Also  since  the 
arguments  are  large  in  most  of  the  region  the  Airy  functions 
can  be  approximated  by 


At 

(’-£’)  OC  Stvt  1 

(Uw+  ?) 

'Em 

(-  £ ) oc  si VI 

By  keeping  track  of  the  phase  differences  of  the  sines  between 
the  two  boundaries  of  each  layer  and  adding  them  up  over  the 
whole  water  depth  one  can  approximately  determine  the  number 


A 
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of  half-waves  accomodated  in  the  water  column  and  thus  de- 
duce the  wave  number  for  the  highest  mode  and  at  the  same 
time  th:e  total  number  of  modes  supported  bv  the  profile: 

M 


or  using  (3*6) 

\/2  to 


A/s  2 


37T  /c  .Va-  *r-t  A: 

|S]  J 


9b) 


Usually  a good  estimate  can  be  obtained  by  approximating 
the  profile  by  two  or  three  straight  lines. 


With  (3.7b)  one  obtains 

f r Ci 


•him- 


-►  £ I CA5+ Ab)V!--  (AjtAfc-  1+  X I | } . 


For  A& 


N S j A* «•  4* »A1» 


‘tw» 


+ «,  At*-  J 


~ H , 


iC  • ^/-1- 

wvn 


j 
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where  H is  the  total  water  denth.  For  the  parameters  chosen 
above  one  obtains  843  modes.  For  = 0 the  result  would  be 

( ^ 3 ) « 

N « M-  ^ A.  [A‘/l- 

3 ( ***  (3.12) 

\ s 215-  •».<<•• . 

Two  approximations  of  (3.12)  are  of  special  interest  : 


A, 


and 


A/s 

4Vi 

f 

3 

<WnVa- 

N =• 

lif£ 

f 

3 

«r  . Va. 

Hence  it  turns  out  that  the  total  number  of  modes  lies 
approximately  in  the  range 


N - ( * * "i  ) 


4/1  f 


-c. 


H VA 


(3.13) 


with  A = maximum  sound  speed  difference , H = total  water 
depth . 

To  determine  the  average  spacing  between  eigenvalues 
one  only  has  to  divide  the  difference 


— •'  - T 
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by  the  total  number  of  eigenvalues  as  determined  above. 
However,  some  care  is  necessary  to  obtain  a useful  result. 
It  turns  out  that  in  the  two  regions 


*>Vci**  < Uz  < , 

(3.15) 

o'/cj-  * kz  < w‘/C  ; 

(3.16) 

the  eigenvalues  are  spaced  very  differently.  For  the  first 

region  to  which  the  eigenvalues  are  confined  when  ct  = c 

b max 

(i.e.  Aj,  = 0)  the  approximate  structure  of  the  eigenvalue 
spectrum  can  be  found  by  calculating  the  number  of  modes  as 
a function  of  the  wave  number  for  a two-layer  profile 
(i.e.  k^  = 0).  Then  the  total  phase  accumulated  in  the  pro- 
file is  given  by  the  two  Airy  function  arguments  at  z = , 

A te)  = A n^). 


If  one  takesA(Zg)  as  a function  of  k varying  between  0 and 
one  can  obtain  the  number  of  modes  as  a function  of  k as 


'Yl  = 


Vz 


jL  5 

■r-  **  I 


L , A 


with£(k)  determined  by 


As 


tod 


to 


+ A ( ki) 
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or 


A(W)S  (<J*/'»)Sk*(eJ./w)(k^ -k). 


(3.17) 


Hence 


OP. 


(3.18) 


or 


V"V> 


JL  - L /,.  ■>  _ / 3 V4  JT-f ' Vj  / t . V 

k-«  ‘T)  (^fAs) 

i.e.  the  wave  number  changes  with  the  two-thirds  power  of 
the  mode  number.  Differentiating  with  respect  to  the  mode 
number  yields 

dt/d-n.  « - | v3  ov" 1/J  Alc/A-*-. 


For An  = 1 one  obtains  the  spacing  between  eigenvalues 
a function  of  the  mode  number: 


as 


OP 


-1/S 


.(3.19) 


For  the  region  (3.16)  the  eigenvalue  distribution  is  quite 
different  as  one  can  see  by  inverting  (3.11)  and  differentiating 
it  with  respect  to  n,  after  replacing  by  (3.17): 


^ H [«-wV->  ( k.,  - k ) ]*'*  ( 3 . 


/C 

with  the  result: 

= - (Tr<mU/lfHx)nl 


20) 


(3.21) 
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i.e.  in  this  region  the  snacing  between  eigenvalues  increases 
linearly  with  mode  number,  whereas  in  the  region  k > 
it  decreases  slowly.  Hence,  the  average  eigenvalue  spacing 
^k>  is  a useful  quantitv  only  in  the  latter  region,  since 
there  the  spacing  remains  relatively  constant  except  for 
the  first  few  modes.  can  be  defined  for  this  case 

by  dividing  (3.14)  with  = 0 bv  (3.12): 


Taking  as  a special  case  =-A.^  and  Ax  , one  gets 


<Ak>» 


-is -(C)" 


(A2  A s ) . 


Hence  a good  approximate  guess  is 


(3.23) 


where  H is  the  total  water  denth. 

This  expression  is  useful  for  establishing  the  step 
size  for  the  wave  number  in  the  search  for  the  eigenvalues. 
Upon  some  further  manipulation  one  can  obtain  the  expression 


A l<  M ^ 2 

/ yi*  \-*/i 


- z y n 

'31a//  . 


£Ak)  “ i U 


(3.24) 
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From  this  one  can  see  that  near  the  highest  mode  (n  = H)  the 
actual  spacing  is  2/3  of  the  average,  for  n = N/8  i.e.  for 
a very  low  mode  the  actual  spacing  is  4/3  of  the  average. 

Hence  by  using  6 x^k>  as  the  step  size  one  has  approximately 
four  steps  between  the  highest  modes , 16  steps  between  modes 
N/64  and  N/64  + 1 which  is  usually  sufficient  to  catch  all 
modes  except  when  nearly  degenerate  modes  are  present. 
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IV.  SOLUTION  OF  THE  CHARACTERISTIC  EQUATION 


As  outlined  in  section  2,  the  eigenvalue  problem  (1.1) 
can  be  recast  into  the  form  (2.29)  or  (2.30)  with  the  pro- 
vision that  the  functions  u”(z,?)  and  ub  (z.f  ) or  u , fz,?  ) 

P P » P*,*  ' 

satisfy  the  wave  equation  (1.1a)  and  the  boundary  conditions 

(1.1b)  and  (1.1c).  This  requirement  has  been  met  by  starting 

solutions  of  (1.1a)  at  the  boundaries  z = 0 and  z 

matching  the  boundary  conditions  there  and  propagating  them 

layer  by  layer  towards  the  joining  point  z = z (for  2.29) 

B 

or  z = zM  (for  (2.30)). 


The  subscripts  p on  the  solutions  indicates  that  the 

particular  functions  describe  the  pth  eigenmode  corresponding 

to  the  eigenvalue  k^.  Hence  all  the  solutions  u”,  ub  and 

can  be  considered  implicit  functions  of  a variable 

k which  for  the  pth  eigenfunction  assumes  the  value  k = k . 

P 


The  task  of  solving  the  eigenvalue  problem  (1.1)  is  thus 
reduced  to  finding  those  values  of  k for  which  (2.29)  or  (2.30) 
are  satisfied,  more  precisely  to  finding  the  zeros  of  the 
functions 


v( 


Pb 


(4.1) 


or 


—j 


S 
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W(M=  u+  (mM  ,flU)  u'  (*„,?,!*) 

(aH/P’/kl-  <4-2) 

As  eq.  (4.2)  suggests,  finding  the  eigenvalues  corresponds 
to  finding  the  zeros  of  the  Wronskian  between  the  down  and  up 
solutions  at  the  joining  point.  Eq.  (4.1)  is  essentially 
equivalent  to  (4.2)  except  that  in  this  case  use  has  been  made 
of  the  fact  that  the  up-solution  is  a simple  exponential. 

This  made  it  possible  to  factor  out  the  up-solution  and  in 
fact  drop  it  since  it  is  nonzero  everywhere. 

An  important  property  of  the  function  W(k)  is  that  it 
alternates  signs  between  successive  eigenvalues.  For  (4.1) 
this  fact  can  be  deduced  in  a straightforward  way  if  this 
equation  is  interpreted  as  a matching  condition  of  the 
derivatives  of  the  down-  and  up-solutions : 

u.wfz£,?,  t)'-  .P,  k)  ' 

where  subsequently  the  relation 

u8(Zi , f , y(?,  k) 

and  the  matching  condition  of  the  functions  themselves 


have  been  used 
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It  is  well  known  that  eigenfunctions  of  the  type 
Up(z,f  ),  regardless  of  their  specific  form,  always  have  a 
standing  wave  type  pattern  of  some  sort  with  p nodes  for 
the  pth  eigenmode.  Hence  if  the  functions  u”  (z,^ ) are 
always  made  to  start  out  with  the  same  slope  at  one  boundary j 
e.g.  the  surface,  they  will  have  alternately  positive  and 
negative  values  at  the  other  boundary.  Therefore,  since  there 
are  no  modes  within  the  ocean  floor,  uw  (zb,  , k)  , and 

v - 

because  of  (2,23)  also  u (zb,  ^ , k)  , have  different  signs 

between  successive  eigenvalues.  Hence  the  derivatives  of 
w 

u (z^,  <$>  , k)  as  a function  of  k oscillates  continuously 
between  positive  and  negative  values  whereas  the  derivatives 
of  ub  Cz^,  £ , k)  is  positive  (negative)  as  long  as  uW  (zb,$>  , k) 
is  negative  (positive)  and  it  shows  the  discontinuous  behavior 
indicated  in  Fig.  2,  switching  signs  whenever  uw  (zfe,  ^ , k) 
changes  sign  (which  is  assumed  to  take  place  at  points  a,  b,  c 
etc.  in  Fig.  2),  The  plus  and  minus  signs  in  Fig.  2 mark  the 
areas  where 

Ufc 

The  numbered  circles  indicate  the  locations  where  W(k)  changes 
sign,  i.e.  the  locations  of  the  eigenvalues. 

L 

Fig.  2 


Ic  l«  - 
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For  eq.  (4.2)  an  argument  similar  to  the  one  just  given  apolies 
if  z„  was  chosen  such  that  the  ud  solution  is  of  the  exponential 
type  throughout.  If  both  the  down-  and  up-functions  have 
oscillatory  behavior  at  z = zw  one  can  simulate  their  behavior 
there  schematically  by  the  functions 


U.|  (2,  k)  - A Sin  — - - 2 

£ 


which  have  matching  logarithmic  derivatives  at  z^  whenever  k is 
an  integer.  Then 


[w(k)} 


M 


AJB  sinkvr 


; 


i.e.  , the  Wronskian  is  zero  for  eigenvalues  and  alternate? sign 
only  there.  This  simple  example  can  also  be  used  to  show  that 
had  one  taken  the  condition  of  matching  the  logarithmic  deriva- 
tives directly  j 

V ( k ) - , 

Uj,  Ut 

instead  of  in  the  form  of  the  Wronskian  (4.2),  the  desirable 
property  of  the  signs  alternately  only  at  the  eigenvalues 
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I 

i 


would  be  lost: 


While  V(k)  still  has  zeros  only  for  k = integer  the  sign  of 
V(k)  now  switches  not  only  when  an  eigenvalue  is  passed  (i.e. 
for  k = integer)  but  also  when  one  of  the  sines  in  the  denomi- 
nator has  a zero. 

Even  though  the  arguments  given  above  were  not  rigorous 
the  qualitative  conclusions  apply  directly  to  the  general  case 
of  arbitrary  eigenfunctions  as  long  as  these  eigenfunctions 
are  all  non-degenerate.  If  this  restriction  is  lifted  and 
degenerate  eigenmodes  are  possible  the  situation  may  become 
rather  complicated  and  the  simple  method  of  finding  the  eigen- 
values outlined  in  this  section  would  certainly  no  longer  work. 
However,  it  is  very  unlikely,  as  long  as  one  is  dealing  with 
empirical  sound  velocity  distributions,  that  one  should  encounter 
cases  with  truly  degenerate  eigenmodes.  On  the  other  hand, 
almost  degenerate  eigenmodes  may  occur  occasionally  in  cases 
where  the  sound  velocity  distribution  has  pronounced  local  minima. 
The  method  of  this  section  may  then  still  be  applied.  However, 
a little  care  may  be  necessary  to  catch  the  sign  change  of  W(k) 
between  the  almost  degenerate  eigenvalues. 
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For  the  purpose  of  this  report  it  will  be  assumed  that 
W(k)  changes  sign  at  every  eigenvalue  and  only  there.  Then, 
in  order  to  find  the  eigenvalues,  i.e.  those  values  of  the  k 
for  which 

r* 

W(k 1=0 


is  satisfied,  one  may  simply  use  the  following  procedure.  W(k) 
is  calculated  as  a function  of  k with  k changing  in.  small  enough 
steps  and,  whenever  the  result  changes  sign  between  two  successive 
steps,  one  switches  over  to  an  iteration  routine  which  locates 
the  exact  position  of  the  zero.  The  range  of  discrete  eigen- 
values is  restricted  to 

' (4.3) 

The  lowest  eigenmode  is  found  for  k*k  , Hence  the  search  is 

max 

started  at  k = kmax  and  proceeds  by  decreasing  k steps  of  XKDj 
the  quantity  XKD  has  to  be  chosen  properly.  If  it  is  too  small, 
excessive  amounts  of  computer  time  will  be  required  to  complete 
the  search.  If  it  is  too  large  the  search  will  locate  only 
part  of  the  eigenvalues.  A good  choice  is  usually 

XKD  & 4 /6 / (4 .4) 

where  ^k^  is  the  average  eigenvalue  spacing  as  defined  in 

section  3.  If  the  bottom  velocity  c is  much  larger  than  c 

o max’ 

the  maximum  sound  velocity  in  the  water,  one  may  save  consider- 
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able  computation  time  by  lettering  XKD  increase  linearly 
with  mode  number  for  wave  numbers 

l<  < k . S Ui/fC 

Occasionally  there  are  almost  degenerate  eigenmodes  which  are 
separated  by  only 

/400 


or  even  less.  In.  those  cases  it  is  usually  still  best  to  use 
(4.4)  and  catch  the  missed  eigenvalues  by  a more  detailed  local 
search. 

Formally  eq,  (4.1)  and  the  method  outlined  for  finding 
u (z,^>  , k)  is  sufficient  to  find  all  eigenvalues.  In  practice, 
however,  several  difficulties  arise  which  have  to  be  dealt  with. 
Usually  no  problems  arise  for  intermediate  or  high  mode  numbers 
since  in  these  cases  the  function  W(k)  is  sufficiently  well- 
behaved  near  the  sea  floor  z = z^»  The  problems  which  do  appear 

are  basically  caused  by  the  fact  that  for  low  mode  numbers 

w ** 

u (z,  £ , k)  has  only  a small  oscillatory  part  near  the  sound 
speed  minimum  and  two  long  exponential  sections  between  there 
and  the  surface  and  the  ocean  floor,  within  which  the  function 
drops  sometimes  over  several  hundred  orders  of  magnitude  (see 
section  3).  However,  despite  the  exponentially  falling  appearance, 
both  exponentially  falling  and  exponentially  rising  components 

hvfp  to  be  present  in  the  solutions  to  allow  matching  of  the 

logarithmic  derivatives  at  the  boundaries.  Very  near  the  eigen- 
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values^  the  coefficients  for  the  exponentially  rising  parts 

are  of  course  extremely  small  so  that  the  required  exponentially 

falling  appearance  can  be  achieved.  However  for  wavenumbers 

away  from  the  eigenvalues  matching  at  the  boundary  between  the 

oscillating  and  the  exponential  part  of  uw  require  large 

coefficients  for  the  exponentially  rising  part.  This  way  W(k) 

can  assume  very  large  values  near  the  sea  floor  which  may  easily 

exceed  the  capacity  of  the  computer.  And  more  seriously,  since 

a zero  of  W(k)  has  to  achieved  essentially  by  cancellation 
w 

between  u an(j  derivative  the  loss  of  accuracy  due  to  very 
near  cancellation  of  very  large  numbers  does  not  permit  to 
obtain  a solution.  Even  the  use  of  double  precision  arithmetic 
is  not  sufficient  in  many  cases.  Four  steps  were  taken  to 
permit  the  calculation  of  the  eigenvalues  despite  these  diffi- 
culties : 

(1)  Instead  of  always  matching  the  down-  and  up-solutions 
at  the  sea  floor  (i.e.  using  eq,  (4.1))  the  matching  is  sometimes 
carried  out  at  some  appropriate  boundary  between  two  layers  in 
the  water  (i.e.  eq.  (4.2)  is  used).  This  makes  W(k)  much  more 
well-behaved  for  the  low  mode  numbers. 

(2)  In  the  cases  which  up-layers  are  introduced  the 
solutions  increase  very  rapidly  between  the  sea  floor  and  the 
matching  point.  To  prevent  exceeding  the  capacity  of  the 
computer  the  function  Uj(z,^>  ,k)  and  its  derivative  u^^2*  f»  , k) 
are  scaled  down  at  each  layer  boundary  in  steps  of  10 10  until 

-V-*  10"5- 
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(3)  If  the  solution  is  of  the  exponential  type  over 
the  whole  extent  of  a wide  layer  the  propagation  of  the 
solution  between  the  two  boundaries  becomes  inaccurate  if 

it  is  carried  out  in  one  step.  In  these  cases  the  layers 
are  divided  into  sublayers  and  the  solutions  are  only  propa- 
gated from  one  sublayer  boundary  to  the  next.  This  way  the 
chances  for  errors  due  to  near  cancellation  of  large  numbers 
is  very  much  reduced. 

(4)  If  the  Airy  function  arguments  are  large  and 
positive  the  corresponding  values  for  the  Airy  functions 
can  easily  exceed  the  capacity  of  the  computer.  In  these 
cases  the  Airy  functions  can  be  approximated  by 

Ai(?  ) ~ *** 

a UVA 

3;  (?■)  «*• 

with  corresponding  expressions  for  the  derivatives,  where  £ 
may  be  of  the  order  of  100  or  even  more. 

Fortunately  it  turns  out  that  one  can  find  a formu- 
lation of  the  problem  (see  eqs.  (2-15)  and  (2-16)  where  these 
very  large  or  very  small  numbers  are  never  needed  individually 
but  only  in  pairs  in  such  a way  that  each  pair  contains  one 
positive  and  one  negative  exponent.  Hence,  by  carrying  fa>, 
g( £ ) etc.  and  the  exponents  along  separately  one  can  avoid 
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overflows  and  underflows  since  the  exponents  appearing  in 
the  results  are  always  only  the  small  differences  of  two 
large  numbers . 

The  computer  code  described  in  this  report  uses  the 
method  for  solving  the  characteristic  equation  outlined 
earlier  in  this  section.  With  the  additional  features  just 
discussed  it  can  be  used  to  solve  a wide  range  of  practical 
problems.  Propagation  problems  in  the  deep  ocean  (H^  4000  m) 
can  be  solved  for  frequencies  up  to  several  hundred  Hz  with 
high  precision.  For  shallow  water  cases  Cwater  depth  = a few 


hundred  meters)  the  frequency  range  in  which  the  program  can 
be  used  extends  up  to  several  kHz.  If  a few  precautionary 
measures  are  taken  in  preparing  the  profile  data  (see  Section  7) 
then  the  range  of  application  is  basically  only  limited  by 
the  amount  of  computer  time  one  is  willing  to  spend. 

The  computer  code  consists  of  the  main  routine  SEARCH 
and  a set  of  subroutines,  A detailed  description  of  these 
programs  is  given  in  the  next  two  sections.  It  should  be  noted 
that  use  has  been  made  of  this  program  in  the  following  publi- 
cation: A.  Nagl,  H.  Uberall , A.  Hang  and  G.  L.  Zarur,  "Adiabatic 
Mode  Theory  of  Underwater  Sound  Propagation  in  a Range  Dependent 
Environment",  J.  Acoust.  Soc.  Am.  63,  73? (1978). 
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V.  DESCRIPTION  OF  THE  COMPUTER  CODE  (MAIN  PROGRAM) 

Purpose  of  this  section  and  the  next  is  to  discuss 
the  set  of  computer  codes  written  to  solve  the  problem  out- 
lined in  section  2 using  the  method  described  in  section  4. 

The  set  consists  of  the  main  program  SEARCH  to  which  this 
section  is  devoted  and  the  subroutines 

PROP,  LARGE,  START,  PR,  LG,  DAIRY,  EXTREM,  SINCOS, 
PARAM  and  SELAUT 

which  will  be  described  in  the  next  section. 

The  main  program  obtains  input  data  from  the  sub- 
routine PARAM,  which  provides  mainly  control  parameters, 
and  the  file  NFILE(A)  in  which  the  sound  velocity  profiles 
are  stored.  The  inputs  consist  of  the  quantities  (indices 
I refer  to  ith  layer  or  ith  layer  boundary,  indices  J refer 
to  jth  profile): 

NUMBER  = number  of  boundaries  for  selected  profile  j 
(X(I),  Y(I)  , 1 = 1,  NUMBER)  = depths  z and  sound 
velocities  cCz^/f  ) for  profile  selected  j 
CBOTM  s sound  velocity  within  ihe  sea  floor.  These  data  come 
from  the  file  NFILE(l),  the  following  ones  from  the  subroutine 


PARAM 
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FREQ  = sound  frequency  f 
RH01  = water  density  ^ 

RH02  = density  of  sea  floor  material  9 D 
XKD  = step  size  of  wave  number  k 
IPRMIN  = first  profile  to  be  investigated 
IPRMAX  = last  profile  to  be  investigated 
XKST(J)  = starting  value  of  k for  eigenvalue  search 
for  jth  profile 

XINCR1,  = depth  increment  and  maximum  depth  at  which 
XMAX1 

depth  functions  are  to  be  calculated 

NOPT  = parameter  to  indicate  whether  or  not  the  depth 

function*  are  to  be  calculated  at  a second 

set  of  depths.  (NOPT  = 1,  NO;  NOPT  = 2,  YES) 

XINCR2 , = depth  increment  and  maximum  deoth  for  second 

XMAX2. 

set  of  depth  function  values  (needed  only 
if  NOPT  =2), 

NFILE(K)  = Names  of  data  files 

NFILE(l)  = Profile  data  file 


NFILE(2)  = Intermediate  file 

NFILE(3)  = Depth  function  values  at  source  depth 
NFILE(4)  = Depth  function  values  at  receiver 
depth 

NFILE(5)  = Eigenvalues 

NFILE(6)  = Wave  functions  and  derivatives  at 


layer  boundaries 
NFILE(7)  = Intermediate  file 
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JMIND  = maximum  number  of  up-layers  allowed  (once  the 

, down-solution  has  reached  this  layer  a switch 

to  the  up-solution  part  may  occur. ) 

ICT  = maximum  number  of  eigenvalues  allowed  for  this  run. 

NPPL(I)  = number  of  depth  function  values  in  ith  layer. 

ISLAY  = number  of  points  at  which  the  depth  functions  are 

calculated  (dimension  of  first  depth  function  array/. 

ISOU  = index  of  element  in  depth  function  array  for  source 
depth. 

IREC1  = index  element  in  depth  function  array  for  receiver 
depth  for  those  modes  for  which  there  are  no  up- 
layers  (if  there  are  up-layers  the  index  is  called 
IREC  and  is  computed  in  the  program) , 

IRLAY  = dimension  of  second  depth  function  array  for  those 
modes  for  which  there  are  up- layers  (needed  only  if 
NOPT  =2). 

IRLAY1  = dimension  of  second  depth  function  array  for  modes 
without  up-layers  (needed  only  if  NOPT  = 2). 

In  executing  the  main  program  a number  of  essential  steps  can 
be  distinguished  which  are  described  below.  In  order  to  facilitate 
associating  variables  used  in  the  previous  section  with  the  corre- 
sponding symbols  in  the  computer  code  both  labels  are.  sometimes 
given,  connected  by  an  arrow.  In  the  program  listing  provided  in 
the  appendix  the  instruction , or  set  of  instructions , each  of  the 
following  steps  is  concerned  with,  are  marked  by  arrows  or 
brackets  and  are  labeled  by  the  appropriate  step  numbers. 
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(1)  Control  parameters  are  read  in  from  PARAM 

(2)  A profile  is  selected  and  the  data  describing  it  are 
read  in. 

(3)  The  wave  number  k ->  XK  for  the  first  iteration  is  deter- 
mined. Generally  the  starting  value  is  . 

However,  if  for  a particular  profile  a value  XKST 

is  provided  the  iterations  will  start  with  this  value. 

(4)  This  step  determines  the  slopes  S(I)  and  the  internal 
derivations  for  each  layer.  All  quantities  are  multi- 
plied with  powers  of  H such  that  they  come  out  dimension- 
less (H  s total  water  depth) . 

(5)  The  search  for  the  first  eigenvalue  is  started  at 
deterrent  number  890.  The  program  returns  to  this 
point  each  time  the  calculations  for  a particular  eigen- 
value have  been  completed  in  order  to  initiate  the  search 
for  the  next  eigenvalue. 

(6)  At  statement  number  990  the  wave  number  k XK  is 
reduced  in  steps  of  Ak  XKD.  For  each  value  of  k thus 
selected  the  program  proceeds  to  calculate  the  Wronskian 
W(k) . This  process  continues  until  W(k)  has  different 
signs  for  two  subsequent  values  of  k.  Then  the  value 

of  k is  determined  by  a different  part  of  the  program  as 
described  below.  Control  of  the  value  of  ,k  remains 
there  until  k has  converged  to  an  eigenvalue. 

(7)  Calculation  of  W(k) , whatever  the  value  of  k,  is  commenced 


at  the  statement  number  994.  First  y-*  UK  is  calculated. 
This  is  up  to  a minus  sign  the  logarithmic  derivative  of 
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the  function  u (z)  inside  the  sea  floor. 

P 

(8)  Then  in  a DO-loop  which  extends  to  statement  3 the 

solution  in  the  water  is  calculated  by  starting  at  the 
surface  with  the  boundary  conditions  u”  (0)  = 0 and 
u”  (b)'  = letl/lfand  propagating  the  function  layer  by 
layer  down  towards  the  sea  floor.  The  solution  and 
its  derivative  at  the  end  of  a layer  (UOUT,  UPOUT) 
in  terms  of  the  solution  and  its  derivative  at  the 
beginning  of  the  layer  (UIN,  UPIN)  are  calculated  in  one 
of  the  three  subroutines  START,  LARGE,  PROP.  PROP  is 
chosen  if  the  argument  £-»DX  of  the  Airy  functions  is 
smaller  or  equal  DM,  which  is  chosen  to  be  11.  LARGE 
is  called  if  the  argument  is  larger  than  DM.  START  is 
used  for  the  first  layer  regardless  of  the  argument. 
However,  START  can  call  LARGE  if  the  argument  is  larger 

than  15  in  at  least  part  of  the  first  layer.  Similarly 

call 

PROP  will /LARGE  if  the  argument  is  ^ 11  at  the  beginning 
of  the  layer  but  grows  > 11  somewhere  within  the  layer. 
Conversely,  if  the  argument  drops  from  ) 11  to  4 11  some- 
where within  a layer,  LARGE  calls  PROP  at  this  point  to 
complete  the  layer. 

This  process  is  carried  through  all  the  way  to  the 
sea  floor  unless  at  the  start  of  the  calculations  for 
some  layer  it  is  found  that  the  argument  of  the  Airy 
functions  at  the  end  of  the  layer  is  larger  than  7 , 

^ and  that  in  addition  the  layer  number  J is  larger  than  JD, 


i 
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the  number  of  completed  down  layers , which  can  be 
set  as  a parameter.  In  this  case  control  of  the 
program  is  transferred  to  another  section  of  the  program 
which  is  described  in  the  next  two  paragraphs. 

If  this  transfer  does  not  occur  and  the  propagation 
of  the  solution  through  all  the  layers  is  completed  the 

n 

program  proceeds  to  calculate  u (zn) ->SCALE  and  W(k) 

p B 

FUNC  which  in  effect  corresponds  to  the  Wronskian  of  the 

' 

sea  floor  as  outlined  in  section  4 , and  then  skips  to 

statement  45  in  order  to  calculate  the  overall  normalization 

constant  (see  paragraph  12.)  The  program  actually  assumes 

Up  (zb}  = 1 and  up  (zB)f  ~~¥  with  UK.  In  this 

sense  SCALE  is  the  mismatch  between  uW  (z.  ) and  uB(z„). 

p b p B 

• w 

Since  u (z)’  is  calculated  as  the  derivative  with  respect 
P 

to  the  argument  of  the  Airy  function  (see  Section  2)  and 

not  with  respect  to  z,  uB  (zD) '-*UPBOTM  has  to  be  divided 

P B 

by  the  internal  derivative  at  z = zg  (=  -AL)  before  it 
can  be  used  to  calculate  the  Wronskian.  The  additional 
factor  is  contained  in  the  definition  of  UPBOTM  to  keep 
the  derivative  dimensionless.  This  is  also  necessary 
since  u^  (z)’  is  dimensionless  by  virtue  of  the  fact  that 
the  argument  of  the  Airy  functions  is  defined  to  be 
dimensionless. 


(9)  The  section  of  the  program  described  in  this  and  the 

next  paragraph  are  entered  if  the  arguments  of  the  Airy  fns 
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grow  too  large  in  the  lower  layers  of  the  profile  for 
a reliable  solution  of  the  eigenvalue  problem  to  be 
• possible  with  the  procedure  outlined  in  the  previous 
paragraph.  If  this  part  of  the  program  is  entered  the 
result  of  the  downward  section  of  the  calculation  (UIN, 
UPIN)  is  kept  for  later  use  and  a new  solution  is 
started  at  z = po  . The  coefficient  of  the  solution 
within  the  sea  floor  is  chosen  such  that  the  value  of 
the  solution  at  the  sea  floor  is 

Up  = XFACT  * 10" 

The  quantity  XFACT  is  at  first  set  equal  to  one.  Later 

it  is  adjusted  to  provide  a match  of  the  upward  part 

of  the  solution  with  downward  part.  Since  u (z)  for 

P 

z ^ Zg  is  an  exponential  the  derivative  at  the  sea  floor 
is  simply: 

uj?  (*“*.)'  =-  UK*  XFACT  . 

In  the  program  this  result  is  divided  by  (-AL)  in  order 
to  give  the  derivative  with  respect  to  the  argument  of 
the  Airy  functions  rather  than  with  respect  to  z.  The 
reason  for  this  was  explained  in  section  2,  The  factor 
H (= total  water  depth)  which  also  appears  is  included 
to  keep  the  derivative  dimensionless. 


UDIN=  RHOZ/UHOl  * XF4CT  * lo-ir 

Starting  with  these  values  the  solution  in  the  water, 
uw(z),  is  then  propagated  upward  layer  by  layer  in  a 
DO-loop  which  extends  to  statement  #4  whereby  the  solu- 
tion and  its  derivative  at  the  end  of  a layer,  i.e.  the 
upper  boundary  of  the  layer  (UDOUT,  UPDOUT)  are  calcu- 
lated in  terms  of  the  solution  and  its  derivative  at  the 
beginning,  i.e.  the  lower  boundary,  of  the  layer  (UDIN, 
UPDIN).  To  carry  out  this  calculation  the  subroutines 
PR  and  LG  are  called  if  the  arguments  of  the  Airy  function 
are  ^ 11  and  ^ 11  respectively.  Transfers  between  PR 
and  LG  within  a layer  can  occur  analogous  to  the  procedure 
outlined  in  paragraph  8. 

Since  for  the  lower  modes  the  upward  solution  may  grow 

3 /2 

very  rapidly  over  many  orders  of  magnitude  (OCexp  2/3  If 
sometimes  % /-w  50  or  even  more)  it  was  found  necessary 

to  provide  the  possibility  of  scaling  down  the  solutions 
and  their  derivatives  at  the  end  of  each  upward  layer. 
Generally,  the  downscaling  consists  of  reducing  u and  u’ 
in  steps  of  10^  until  |u  [ < 10  \ No  scaling  is  performed 
in  the  uppermost  up-layer.  To  keep  track  of  the  number 
of  scalings  at  the  end  of  each  layer,  this  number  is 
stored  in  the  array  ISC(J).  This  array  is  updated  for 
each  new  value  of  k.  However,  once  the  Wronskian  has 
changed  signs  the  number  of  scalings  at  each  boundary 
is  kept  constant  in  order  to  prevent  the  scaling  pro- 
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cedure  from  disrupting  the  iteration  process  which  finds 
the  k for  which 

V (IO  = o . 

If  the  number  of  scalings  in  two  successive  iterations 
were  allowed  to  be  different  the  convergence  of  the 
iteration  process  would  be  uncertain. 

The  fact  that  the  quantity  XFACT  may  change  its  value 
adds  another  complication.  XFACT  corresponds  to  the 
mismatch  between  the  down-and  up-solutions  at  their 
joining  point.  Until  k has  converged  to  some  eigenvalue 
kp,  XfACTis  actually  kept  equal  to  one.  This  does  not 
interfere  with  finding  the  eigenvalue.  However,  in  order 
to  get  a continuous  solution  and  to  find  the  correct  value 
for  the  normalization  integral  it  is  necessary  at  the  end 
to  recalculate  the  solution  with  the  eigenvalue  just  found 
and  with  XFACT  assuming  its  proper  value. 

In  this  calculation  the  number  of  scalings  in  each  layer 
may  be  different  than  before.  Hence  in  order  to  keep 
track  of  the  number  of  scalings  in  the  final  iteration 
another  array  has  to  be  introduced:  MSC(J).  If  the 
solutions  are  required  at  two  different  set  of  depths 
(NOPT  = 2)  an  additional  iteration  is  carried  out  where 
the  number  of  scalings  to  be  performed  in  each  layer  is 
taken  from  MSC( J) . 

-J 
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(11)  After  the  calculation  for  all  the  up-layers  has  been 

completed  the  Wronskian  at  the  layer  boundary  where  the 
up-  and  down-solutions  join,  is  determined.  If  the 
program  is  in  the  phase  where  an  eigenvalue  has  been  found 
in  the  previous  iteration  (MFLAG  = 1)  two  additional  cal- 
culations are  carried  out  in  this  section.  First,  the 
array  NSCCICOUNT,  K)  is  determined  which  gives  the  total 
number  of  scalings  which  have  been  carried  out  on  the 
solution  with  mode  number  ICOUNT  between  the  sea  floor 
and  the  Kth  layer.  In  other  words,  in  order  to  obtain  the 
correct  values  for  the  solution  for  the  mode  p -*  ICOUNT 
in  the  Jth  layer,  the  solution  as  stored  in  the  file  at 
this  point  has  to  be  multiplied  by 
10*  N SO  ( ICOUNT t K) 

The  second  operation  carried  out  in  this  section  concerns 
the  contributions  of  the  up-layer  to  the  normalization 
integral.  They  are  also  affected  by  the  scalings  and 
have  to  be  corrected.  XS(J)  and  XT(J)  are  contributions 
to  the  normalization  integral  accumulated  in  layer  J while 
the  solution  was  calculated  in  subroutines  PR  and  LG, 
respectively.  It  was  assumed  that  only  the  uppermost 
up-layer  (J=JX)  contributes  significantly  to  the  normali- 
zation integral.  The  solution  can  be  assumed  to  be 
exponentially  decaying  for  at  least  a considerable  part 
of  this  layer  so  that  in  the  lower  layers  the  square 
of  the  solution  will  be 
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insignificantly  small.  Hence,  the  contribution  of  all 
the  up-layers  to  the  normalization  integral  was  taken  to 


XNORM  2 **  ( XS  (XX  ) + XT  (TX))*  10 


-20  *kfS 


where  NS  is  the  difference 


ISC  (T)  - M SC  IT) 


assumed  over  all  up-layers. 

C12)  The  contributions  of  the  down-layers  to  the  normalization 
integral  are  accumulated  through  the  variable  XNORM,  each 
layer  adding  the  quantity 


J I I*"  <*2:  . 


The  contributions  of  the  up-layer  are  calculated  similarly, 
added  up  to  yield  XNORM  as  outlined  in  the  previous  para- 
graph. The  calculation  of  these  integrals  is  described 
in  the  section  on  the  subroutines.  The  contribution  of 
the  solution  inside  the  sea  floor 
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is  called  XN0R1.  It  is  obtained  by  observing  that  up  (z) 
is  an  exponential  with  the  exponent  -ar*  -XK.  Hence  the 


overall  normalization  constant  is  obtained  as 

SaW  = ( + yW0RM  + XWORM  2.)  ^ 

where  the  fact  was  taken  into  consideration  that  SCALE 

w B 

is  the  mismatch  between  ^^u  (zR)  and  ^>Ru^  ^zb^* 
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stored  in  the  file  NTYPE(2)  with  FORMAT  (D). 

(13)  If  MFLAG  = 1 (=  eigenvalue  found  in  previous  iteration) 

and  NOPT  = 2 ( = second  set  of  depth  functions  requested) 

the  program  will  return  to  statement  #994  (see  paragraph  7) 
to  calculate  the  depth  functions  in  the  range 

Oiii  XMaX  2. 

in  steps  of  XINCR2.  These  values  are  also  stored  in  file 
NTYPE(2).  The  format  in  this  case  is  FORMAT  (D  20.8). 

(14)  If  MFLAG  = 1 then  all  calculations  concerning  a particular 

eigenvalue  are  finished  when  statement  #892  is  reached  and 
the  program  returns  from  this  point  to  statement  number 
890  (see  paragraph  5)  to  start  the  search  for  next 
eigenvalue. 

(15)  The  program  enters  this  section  only  once  for  each  eigen- 
value, namely  if  it  is  found  that  the  wave  numbers  for 
two  successive  iterations  differ  by  less  than  10“^.  This 
condition  is  used  as  the  criterion  that  the  search  for  the 
eigenvalue  has  converged.  The  value  k XK  obtained  is 
stored  in  file  NTYPE(2)  with  FORMAT(D).  If  there  are  up- 
layers  the  mismatch  between  the  down-  and  the  up-solutions 

(XFACT)  is  calculated  here  and  the  program  returns  to 
statement  number  994  (see  paragraph  7)  with  MFLAG  = 1 
for  one  (NOPT  * 1)  or  two  (NOPT  =2)  additional  iterations 
in  which  the  following  quantities  are  calculated  and 
stored  (as  already  partially  indicated  above): 

(i)  contributions  from  each  layer  to  the  normalization 


L 


integral 
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(ii)  the  values  of  the  solution  uW(z)  at  one  or  two 

P 

different  sets  of  depths. 

(iii)  the  values  of  the  solutions  and  their  derivatives 

at  each  layer  boundary. 

(16)  This  section  of  the  program  determines  the  wave  number 
k for  the  next  iteration.  If  the  Wronskian  W(k)  has 
not  changed  sign  since  the  last  eigenvalue  was  found 
(IFLAG  = 0)  the  program  returns  to  statement  #990  where 
k-^XK  is  reduced  by  A k — ^XKD . If  the  Wronskian  has 
changed  sign  between  this  iteration  and  the  previous  one, 
JFLAG  is  set  equal  to  1 and  the  wave  number  for  the  next 
iteration  is  calculated  by  quadratic  interpolation: 

If  JFLAG  = 1 already,  k for  the  next  iteration  is 
calculated  by  linear  interpolation: 

W(kJ 

and  the  program  also  returns  to  statement  number  994. 

(17)  This  section  of  the  program  is  entered  after  all  eigen- 
values for  a particular  profile  have  been  found.  It 
essentially  reads  back  all  the  data  written  into  file 


k.*  = k 
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NFILE (1),  pickes  out  the  depth  function  values  at  the 
source  and  receiver  depths  and  normalizes  them.  Then 
it  stores  the  values  at  the  source  depth  for  all  modes 
into  file  FILE  (3),  the  corresponding  values  at  the 
receiver  depths  into  file  NFILE  (4),  the  eigenvalues  into 
fileNFILE  (5),  the  wave  functions  and  their  derivatives 
at  the  larger  boundaries  into  file  NFILE  (6) . In  all 
cases  FORMAT  (86D)  is  used. 


VI.  DESCRIPTION  OF  THE  COMPUTER  CODE  (SUBROUTINES). 

Subroutine  SELAUT 

This  routine  transfers  the  data  for  the  selected  profile 
from  the  profile  data  file  NFILE  (1)  to  the  temporary 
file  NFILE  (2)  converting  them  in  the  process  into  a 
format  which  is  convenient  for  reading  them  to  the  main 
program  SEARCH.  This  arrangement  is  convenient  since 
the  main  program  does  not  have  to  be  changed  if  for  some 
reason  the  arrangement  of  the  data  in  the  profile  data 
field  should  be  altered. 

Subroutine  PARAM 

Purpose  of  this  program  is  to  provide  a convenient  input 
channel  for  all  the  input  parameters  except  the  profile 
data.  Through  the  set  of  input  parameters  provided  it  is 
usually  possible  to  set  up  the  program  for  specific  cases 
without  having  to  change  and  recompile  the  main  program 
and  any  of  the  other  subroutines  every  time.  The  set  of 
parameters  provided  is  listed  in  the  previous  section, 
where  the  main  program  is  described.  Information  on  how 


to  set  up  a run  is  given  in  section  7. 
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Subroutines  DAIRY 

It  calculates  the  Airy  functions  and  their  derivatives 
to  14  significant  digits  for  arguments  in  the  range 
-1000<  DX<  15  . 

It  is  called  from  START,  PROP  and  PR.  The  arguments 
for  which  the  Airy  functions  are  needed  in  these  subroutines 
are  usually  smaller  than  DM,  which  is  11  in  PROP  and  PR 
and  15  in  START. 

If  the  argument  is  larger  than  that  value  in  some  part 
of  a layer  the  program  switches  over  to  LARGE  or  LG. 

The  only  time  when  the  argument  can  be  larger  than  11 
in  PROP  and  PR  is  in  the  case  of  a layer  in  which  the 
argument  decreases  from  above  11  to  below  11.  Then  the 
program  starts  the  layer  in  LARGER  or  LG  and  then  switches 
to  PROP  or  PR  at  tfte  beginning  of  the  sublayer  in  which 
the  argument  decreases  below  11.  At  this  point  the  argu- 
ment may  be  larger  than. 11  by  the  amount  DA  by  which  the 
argument  may  change  within  a sub-layer.  Hence  if  one 
makes  sure  that  DA  ^4  the  arguments  for  which  DAIRY  is 
used  always  stay  below  the  upper  limit  for  which  DAIRY 
is  defined. 

— — S 


For  layers  with  very  small  slopes  in  the  velocity  profile 
the  arguments  may  fall  considerably  below  -1000  (see 
section  3) . For  these  cases  the 

Subroutine  SINCOS 

is  called  which  uses  the  asymptotic  expansion  for  large 
negative  arguments  given  in  (3.1). 

Subroutine  EXTREM 

Calculates  the  Airy  functions  and  their  derivatives 
for  large  positive  arguments.  This  routine  is  called 
from  LARGE  and  LG.  The  arguments  for  which  the  Airy 
functions  are  needed  there  are  usually  larger  than  11. 

An  exception  are  the  calculations  in  a sublayer  in  which 
the  arguments  rise  above  or  fall  below  11,  i.e.the 
sublayers  at  the  beginning  (the  end)  of  which  the  switch 
from  (to)  the  routines  PROP,  START  or  LG  occurs.  Then 
the  arguments  can  be  below  11  by  the  amount  by  which  they 
can  change  within  a sublayer. 

This  means  that  the  ranges  of  arguments  fcr  which  EXTREM 
is  called  extends  from  about  7 to  possibly  several  hundred. 
To  avoid  overflows  resulting  from  the  corresponding  large 
values  for  the  Airy  functions  the  results  are  given  in 
2 parts,  one  being  the  exponent  of  the  asymptotic  form 
(=  ZETA) , the  other  one  the  value  of  the  polynomial 


multiplying  the  exponential  (see  (3.1)).  in  order  to 
take  advantage  of  the  .high  precision  of  the  subroutine 
DAIRY  for  the  smaller  arguments,  this  routine  is  called 
by  EXTREM  to  calculate  the  Airy  functions  for  arguments 
6 15. 

Purpose  of  the  remaining  five  subroutines  (PROP,  PR, 

START,  LARGE,  LG)  is  to  carry  out  the  following  functions: 

(1)  To  provide  the  depth  functions  and  derivative* at  the 
end  of  an  interval  in  terms  of  the  depth  function  and  its 
derivative  at  the  beginning  of  the  interval.  START,  PBOP 
and  LARGE  are  used  for  down-layers,  and  PR  and  LG  for  up- 
layers.  (2)  If  MFLAG  = 1,  to  calculate  the  integral  of  the 
square  of  the  depth  function  over  the  range  of  the  layer. 
(3)  If  MFLAG  =1,  to  calculate  the  depth  function  at 
integer  multiples  of  XINCR.  The  inputs  to  all  5 of  these 
routines  are  with  a few  exceptions: 

AL,  BE  quantities  in  -terms  of  which  the  arguments  of 
the  Airy  functions  are  defined.  According  to 
(2.10)  the  arguments  can  be  written  as 

ff)]  or  - y4L  * (B£>  *), 

where  AL  corresponds  to  the  internal  derivative 
and  BE  to  tip ) j wi'fch.  d ^ ( p)  being 

defined  by  (2.9).  AL  and  BE  can  always  be  assumed 


to  be  constant  within  the  range  (XxXJ)  of  the 
interval.  Hence  the  arguments  of  the  Airy 
functions  are  always  either  a monotonically 
increasing  or  a monotonically  decreasing  function 
of  the  depth. 

X,  XJ  Depth  variables  defining  the  beginning  and  the 
end  of  the  depth  interval  for  which  the  calcu- 
lations are  performed. 

U,  UP  Depth  function  and  its  derivative  at  beginning 
of  the  interval. 

XNEW  Depth  variable  which  is  changed  in  increments 
of  XI NCR  inside  the  subroutines . It  specifies 
the  depth  at  which  depth  functions  may  be  calcu- 
lated and  it  is  used  to  define  the  boundaries 
of  sublayers.  The  input  value  of  this  variable 
is  always  ^ X for  START,  PROP,  LARGE  and  ^ X for 
PR,  LG. 

XINCR  Step  size  with  which  the  depth  variable  XNEW  may 
be  changed. 

TT  Ratio  of  internal  derivatives  at  boundary  between 

two  intervals  as  defined  in  (2.13')  (internal 
derivative  for  previous  interval  divided  by 
internal  derivative  for  this  interval). 


59 


MFLAG  Flag  that  is  set  in  the  main  program  SEARCH  and 
transmitted  to  the  subroutines  through  a COMMON 
statement.  It  indicates  that  an  eigenvalue  has 
just  been  found.  In  the  subroutines  it  sets  up 
the  conditions  for  calculating  the  contribution 
to  the  normalization  integral  from  the  layer  under 
consideration  and,  if  in  addition  MFLAG 1 = 1,  for 
calculating  the  depth  functions  at  prescribed  depths. 

MFLAGCt  This  flag  is  also  transmitted  to  the  subroutines 

through  a COMMON  statement.  It  is  set  to  1 in  the 
main  program  SEARCH  if"  the  maximum  depths  (XMAXl 
or  XMAX2.);  down  to  which  the  values  of  the  depth 
functions  are  requested^ are  within  or  below  the 
layer  under  consideration.  If  the  flag  has  the 
value  1 it  causes  the  subroutines  to  calculate  the 
depth  functions  for  all  values  of  XNEW  which  fall 

inside  the  interval,  provided  MFLAG  also  has  the 

• \ 

value  1. 


t>- * 


60 


The  output  from  the  subroutines  consists  of: 

U^f  l/PJ  Depth  function  and  its  derivative  at  end  of  interval. 

During  execution  of  the  subroutines  this  variable 
has  been  increased  (START,  PROP,  LARGE)  or  decreased 
(PR,  LG)  by  an  integer  multiple  of  XINCR  from  its 
input  value.  Its  output  value  is  always  ^ XJ  for 
START,  PROP,  LARGE,  and  ^ XJ  for  PR  and  LG. 


XNORM  This  is  the  variable  through  which  the  contribution 
to  the  normalization  integral  of  all  the  down-layers 
is  accumulated,  i.e.  whenever  any  of  the  subroutines 
START,  PROP,  or  LARGE  is  called  and  the  condition 
MFLAG  = 1 exists  the  integral  of  the  square  of  the 
depth  function  within  the  layer  boundaries  is 
calculated  and  the  result  added  to  XNORM.  The 
variable  is  transmitted  to  and  from  the  main  program 
via  a COMMON  statement. 

X1SQ,  These  are  the  contributions  to  the  normalization 
X2SQ  integral  calculated  in  the  subroutines  PR  (aXlSQ) 
and  PG  (oC2SQ) . PR  and  LG  have  both  quantities  as 


arguments  since  a transfer  between  the  2 subroutines 
can  occur  within  a layer.  The  contribution  of  all 
the  up-fayers  are  added  up  in  the  main  program. 
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Subroutines  PROP  and  PR  have  a similar  structure  and  are  therefore 
discussed  together.  They  consist  essentially  of  3 sections: 

Section  (1) : Calculations  at  the  beginning  of  the  layer.  The 
coefficients  AJ,  BJ  of  the  Airy  functions  defining  the  depth 
function  in  this  layer  are  determined  here  following  eqs.  (2-13'). 

If  MFLAG  = 1 the  value  of  the  normalization  integral  at  the  start 
of  the  integration  interval  ( =BEG)  is  also  calculated  at  this  point. 

Section  (2):  Calculations  at  sublayer  boundaries.  This  section 
is  entered  in  two  cases: 

(a)  If  the  argument  of  the  Airy  function  at  XJ,  the  end  of  the 
layer,  is  > 11. 

(b)  If  MFLAG i = 1,  i.e.  if  the  calculation  of  the  depth  function 
at  points  within  this  layer  is  requested.  In  both  cases  sublayer 
boundaries  are  established  within  the  layer  at  regular  intervals, 
starting  with  the  input  value  of  XNEW  increased  for  PROP  (decreased 
for  PR)  by  XINCR  and  then  proceeding  in  steps  of  XINCR.  At  each 
boundary  the  values  of  depth  function  (WF)  and  its  derivative  (WP) 
are  calculated.  In  case  (2)  the  depth  function  at  each  sublayer 
boundary  is  stored  in  file  NFILE  (2)  with  FORMAT  (D20.8) . In 
case  (1)  the  program  proceeds  from  one  sublayer  to  the  next  until 
at  some  sublayer  boundary  the  argument  of  the  Airy  function  exceeds 
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11.  Then  subroutine  LARGE  is  called  using  as  inputs 
X = XNEW 
U = WF 
UP  = WP 

TT  = 1 (since  the  internal  derivatives  do  not  change  across 
the  sublayer  boundaries), 

completing  there  the  remainder  of  the  calculation  for  this  layer. 
After  completion  of  the  CALL  statement  the  program  immediately 
returns  to  the  main  program,  if  MFLAG  = 1 the  value  of  the  norma- 
lization integral  at  the  latest  XNEW  is  calculated  (-CEN)  and 
the  quantity  (-CEN  -BEG)  is  added  to  XNORM  before  returning  to 
the  main  program,  i.e.  the  contribution  to  the  normalization  integral 
of  the  region  between  the  beginning  of  the  layer  and  the  XNEW  at 
which  the  transfer  to  LARGE  occurs  as  added  to  the  total. 

Section  (3);  Unless  a transfer  to  subroutines  LARGE  occurred 
in  Section  2 the  calculations  in  PROP  and  PR  are  wound  up  by 
calculating  the  depth  function  and  its  derivative  at  the  end  of 
the  layer  (i.e.  UT,  UPJ)  and  if  MFLAG  = 1 by  calculating  the 
normalization  integral  at  the  end  of  the  layer  (-CEN)  and  adding 
the  difference  to  XNORM.  If  Section  2 was  skipped  (i.e.  if  the 
argument  of  the  Airy  functions  was  smaller  than  11  for  the  whole 
layer  and  if  MFLAG 1 was  0,  so  that  the  depth  functions  did  not 
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have  to  be  calculated  inside  this  layer) , then  XNEW  still  has 
its  initial  value  and  it  has  to  be  updated.  If  the  argument  of 
the  Airy  functions  is  ^ 11  at  the  beginning  of  the  layer  but^>  11 
at  the  end,  only  the  solution  for  the  first  part  of  the  layer  is 
calculated  in  PROP  ( or  PR)  up  to  the  last  sublayer  boundary  at 
which  the  argument  is  less  than  11.  The  solution  for  the  rest 
of  the  layer  is  then  calculated  in  LARGE  (or  LG)  as  indicated  above. 
The  argument  at  which  the  transfer  occurs  may  be  from  very  slightly 
less  than  11  to  considerably  less  than  11.  The  latter  is  the  case 
if  the  slope  is  large  and  the  argument  at  the  next  sublayer  boundary 
is  just  above  11.  The  reverse  case  where  the  argument  is^  11  at 
the  beginning  of  the  layer  and^  11  at  the  end,  i.e.  when  the  calcu- 
lations for  the  first  part  of  the  layer  are  made  in  LARGE  (or  LG) 
and  then  the  transfer  is  made  to  PROP  (or  PR)  is  discussed  in  the 
section  on  LARGE  and  LG. 

Subroutine  START  is  called  to  perform  the  calculations  in  the 
layer  nearest  the  surface.  It  contains  4 sections:  the  first 
three  are  used  if  the  argument  of  the  Airy  functions  remain  below 
15  for  the  whole  layer,  the  fourth  one  if  the  argument  exceeds  15 
for  at  least  part  of  the  layer. 

Section  1:  calculations  at  the  surface.  Calculation  of  the 
coefficients  to  the  Airy  functions  according  to  assumptions 
(2-18)  and  (2-19) . If  MFLAG  = 1 the  value  of  the  normalization 


integral  at  the  surface  is  also  determined  here. 


Section  2:  calculations  at  the  end  of  the  layer.  This 
section  determines  the  depth  function  and  its  derivative 
at  the  end  of  the  layer  and,  if  MFLAG  = 1,  the  contribution 
ZNORM  of  the  layer  to  the  normalization  integral  XNORM. 

Section  3;  updates  XNEW  to  the  largest  value  ^ XJ  and  for 
MFLAGl  = 1 calculates  the  depth  function  inside  the  first 
layer  at  intervals  XINCR  and  writes  them  into  file  NFILE  (2). 
Section  4;  This  section  is  used  to  transfer  the  calculations 
for  the  entire  first  layer  to  LARGE.  This  is  done  if  anywhere 
inside  the  layer  the  argument  of  the  Airy  function  exceeds 
15.  For  the  part  of  the  layer  for  which  the  argument  is 
^15  the  calculation  should  actually  be  done  in  START  but 
it  is  assumed  that  if  the  argument  is  > 15  somewhere  in 
the  layer  it  is  not  much  less  anywhere  else.  This  is  a 
reasonable  assumption  since  surface  layers  with  large  argu- 
ments are  never  very  wide.  So  the  argument  may  be  large  but 
it  cannot  vary  much  over  the  extent  of  the  layer. 

Subroutines  LARGE  and  LG  again  have  similar  structure  and  may 
therefore  be  discussed  together.  These  subroutines  are  used  for 
those  intervals  in  which  the  Airy  function  arguments  are  large 
( y DM) . Such  an  interval  may  extend  over  a whole  layer  or  over 


part  of  a layer.  Since  the  Airy  function  arguments  are  monotonically 

in 

increasing  or  decreasing  with/a  layer,  it  is  clear  that  of  the 
arguments  at  the  beginning  (=  DX)  and  at  the  end  (=  DXl)  at  least 
one  must  be  larger  than  DM,  if  the  routines  LARGE  and  LG  are  involved 
in  the  calculation  for  a particular  layer,  and  one  can  distinguish 
3 cases: 

(1)  DX  > DM  and  DXl  ■>  DM,  The  layer  is  calculated  entirely 
by  LARGE  or  LG  (interval  extends  between  the  layer 
boundaries  X and  XJ.) 

(2)  DX  > DM  but  DXl  <,  DM.  First  part  of  layer  is  calculated 
in  LARGE  or  LG  (interval  extends  from  X to  value  of  XNEW 
at  switch.  For  calculations  of  remainder  of  layer  PROP 
or  PR  are  called. 

(3)  DX  ^ DM,  DXl > DM.  First  part  of  layer  is  calculated 
in  PROP  or  PR.  For  calculations  of  remainder  of  layer 
LARGE  or  LG  are  called  by  PROP  or  PR  (interval  extends 
from  value  of  XNEW  at  switch  to  XJ) . 

All  intervals  for  which  LARGE  and  LG  are  used  are  divided  up 
into  sublayers  of  width  XINCR,  and  the  depth  functions  are 
propagated  from  one  sublayer  boundary  to  the  next  using  eqs.  (2-15) . 
The  first  sublayer  boundary  extends  from  X (=  interval  boundary, 
which  is  the  layer  boundary  in  cases  (1)  and  (2)  and  the  value  of 
XNEW  at  the  switch  in  case  (3))  to  XNEW. 


+ XINCR,  the  first 
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sublayer  boundary  within  the  interval.  Subsequent  sublayers  are 
established  by  changing  the  depth  variable  XNEW  in  steps  of  XINCR 
and  using  these  values  as  boundaries.  The  last  sublayer  extends 
from  the  highest  (lowest)  value  for  XNEW  which  is  4 XJ  ( >-  XJ) , 
or  in  case  (3)  to  the  first  value  of  XNEW  for  which  the  Airy  function 
argument  is  <11. 

Both  LARGE  and  LG  essentially  consist  of  3 sections: 

Section  (1);  At  the  beginning  of  the  interval  only  the  Airy 
functions  for  the  argument  at  this  point  are  calculated. 

Section  (2):  Calculations  for  each  sublayer.  The  solutions  are 
propagated  from  one  sublayer  to  the  next  by  expressing  the  depth 
functions  and  their  derivatives  at  the  end  of  each  sublayer  (Ul, 

UlP)  in  terms  of  the  corresponding  values  at  the  end  of  the  previous 
sublayer  (U,  UP)  and  the  coefficients  following  eqs.  (2-16). 

If  MFLAG  = 1 the  depth  function  value  is  written  into  the  file 
NFILE  (2)  at  every  value  of  XNEW  inside  the  interval. 

If  MFLAG  = 1 the  contribution  of  the  sublayer  to  the  normalization 
integral  is  calculated. 

If  at  some  sublayer  boundary  it  is  found  that  within  the  following 
sublayer  the  Airy  function  argument  goes  below  DM,  the  calculations 
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are  transferred  to  PROP  or  PR,  in  the  case  MFLAG  = 1 after  the 
contribution  of  all  the  sublayers  up  to  this  point  to  the  norma- 
lization integral  is  calculated  and  added  to  XNORM. 

The  argument  at  which  the  transfer  occurs  may  vary  from  slightly 
larger  than  11  to  considerably  larger  than  11.  The  latter  case 
may  happen  if  the  argument  changes  very  much  over  the  extent  of 
the  sublayer  and  it  is  only  slightly  less  than  11  at  the  next 
sublayer  boundary.  This  means  that  the  Airy  function  argument 
could  be  larger  than  11.  This  means  that  the  Airy  function  argument 
could  be  larger  than  15  for  calculations  in  PROP  or  PR.  This  could 
cause  problems  since  the  subroutine  DAIRY  which  is  used  in  PROP 
and  PR  for  calculating  the  Airy  functions  is  accurate  only  for 
arguments  below  15,  and  more  seriously,  the  value  of  the  Airy 
functions  may  exceed  the  capacity  of  the  computer  if  the  argument 
is  allowed  to  go  too  high.  In  order  to  avoid  these  problems  the 
width  of  the  sublayers  should  be  chosen  small  enough  to  keep  the 
change  in  the  Airy  function  arguments  within  a sublayer  below 
4. 

According  to  (3.6)  the  difference  in  the  Airy  function  arguments 
can  be  written  as 

= " ^ [Afen,J - AfO] 

where 


- A (2j  = (A f 4.-/4,)xJA/CR 
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i.e.  the  difference  in  the  A' s is  just  the  sound  velocity  slope 
in  the  layer  times  the  thickness  of  the  sublayer.  Hence 

ArsL  = a//5  (Aif3  AIVCR  . 

So  in  order  to  keep  the  Airy  function  arguments  from  exceeding 
15  for  calculations  in  PROP  and  PR,  XINCR  should  be  chosen  such 
that 


CO 


2/3 


Ai  1 1/1 

17 


Xl*/CR  < if 


(6.4.) 


for  all  layers. 

Section  (3 ) : Calculations  at  the  end  of  a layer  consist  of  propa- 
gating the  solution  from  the  last  sublayer  boundary  to  the  end  of 
the  layer  and  if  MFLAG  = .1  of  calculating  the  contribution  of  the 
vhole  interval  to  the  normalization  integral  and  adding  it  to 
XNORM. 


\ 
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VII.  INSTRUCTIONS  FOR  USING  THE  COMPUTER  CODE 

This  section  outlines  the  steps  necessary  to  obtain  the  solutions 
to  the  eigenvalue  problem  (1.1)  using  the  computer  code  described 
in  the  previous  two  sections. 

(1)  Preparation  of  input  parameters  (entered  into  PARAM) . 

(NFILE (J) , J = 1.5)  Naming  of  files: 

J = 1 profile  file 

J = 2 temporary  file 

J = 3 depth  functions  ai  source  depth 

J = 4 depth  functions  a*t  receiver  depth 

J = 5 eigenvalue  file 

J = 6 file  for  depth  functions  and 

derivatives  at  layer  boundaries 


FREQ,  RHOl,  RH02 
IPRMIN,  IPRMAX 
XKST(J) 


Sound  frequency,  water  density,  seafloor 

density 

First  and  last  profile  to  be  investigated 
Starting  value  of  k for  eigenvalue 
search  for  jth  profile.  Ordinarily 
the  values  of  this  array  are  all  set 
equal  to  1000.  This  way  the  first 
mode  to  be  found  will  be  the  one  with 
mode  number  one.  If  desired  a run 
can  be  started  with  a higher  mode 
number  by  specifying  the  appropriate 
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ICT 


JMIND 


XKD 


values  of  the  wave  numbers  XKST  for 
each  of  the  selected  profiles. 

This  parameter  specifies  the  number 
of  modes  to  be  calculated.  If  a number 
>N  (N  = maximum  number  of  modes  for  a 
profile)  is  selected  all  modes  are 
calculated. 

Maximum  number  of  up-layers  allowed. 
Usually  a good  choice  is  to  take  one 
less  than  the  number  of  layers  below 
the  location  of  the  sound  speed  minimum. 
Step  sifee  of  wave  number  during  search 
for  sign  change  of  W(k).  This  quantity 
determines  to  a large  extent  the  amount 
of  computer  time  required  for  a parti- 
cular profile.  A good  value  to  start 
with  is  XKD  = ^ Ak^ /6  with^k^  as 
defined  by  (3-23) . For  the  region 


k < 


CO 

c 


vc 


max 


» 


c — c 

max  min 


max 

the  spacing  between  eigenvalues  increases 
linearly  with  mode  number  according  to  (3.21) 
Hence,  if 


C 


B 


C 

max 


» 


C 

max 


C . 
min 


with, e.g. , 


C . %.  1490  m/s,  C *2  1540  m/s, 
mm  max 

Cb2£  1800  m/s. 
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XI NCR 1, 
XINCR2, 
NOPT 


[ 


XMAX1 

XMAX2 


r) 


the  spacing  between  eigenvalues  can 

grow  very  large  and  one  can  save  a 

considerable  amount  of  computer  time 

if  the  step  size  XKD  is  made  to  increase 

linearly  with  mode  number  when  searching 

. OJ 

for  eigenvalues  in  the  region  k<“ * 

max 


This  can  be  done  by  making  use  of 
increasing  XKD  by 

TC  c . 


mm 


2 f nz 


(3-21)  . 


\ 


each  time  an  eigenvalue  has  been  found. 

These  quantities  specify  the  depth 
at  which  the  depth  functions  are  calculated. 
In  addition  XINCRl  and  XINCR2  determine 
the  widths  of  the  sublayers. 

If  the  source  and  receiver  depths  are 
related  by  an  integer  or  if  they  have 
a convenient  common  factor  the  depth 
functions  at  both  locations  can  be  calcu- 
lated with  one  pass  through  the  program 
after  a particular  eigenvalue  has  been 
found.  In  this  case  one  sets 
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NOPT  = 1 

and  XINCR2,  XMAX2  are  not  needed. 

XINCR1  is  chosen  such  that  both  the  source 
and  receiver  depth  are  an  integer  multiple 
of  it  and  XMAXl  is  taken  to  be  the  source 
or  receiver  depth,  whichever  is  larger 
(e.g.,  source  depth  = 200m,  receiver  depth 
= 500m,  XINCR1  = 100m,  XMAXl  = 500m) . 

If  the  source  and  receiver  depth  have  no 
convenient  common  factor  (e.g.  sorce  depth 
= 29m,  receiver  depth  = 100m)  the  depth 
functions  at  these  depths  are  calculated 
at  two  different  passes  through  the  program. 
In  this  case 
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NPPL(I),  ISLAY, 
ISOU,  IREC1, 
IRLAY,  IRLAY1 


NOPT  = 2 

XINCR1=  XMAX1  = Source  depth 

XINCR2  = XMAX2  = Receiver  depth. 

The  option  NOPT  = 1 must  also  be  used  if 
the  receiver  position  falls  into  an 
up-layer.  This  case  occurs  when  the 
water  depth  is  very  large  and  the  receiver 
is  located  near  the  ocean  floor.  At 
the  end  of  section  6 it  was  mentioned 
that  there  is  a constraint  on  XINCR, 
which  took  the  form  (6.1).  Accordingly 
XINCRl  and  XINCR2  have  to  satisfy  the 
condition 

XINCRl,  XINCR2  < 1400  f“2/3(  A/h)"1^3 

max 

where  *f  is  the  frequency  and  (A  /h) 

max 

the  maximum  sound  speed  gradient  of  profile, 
are  the  parameters  needed  to  extract 
the  desired  output  data  (eigenvalues, 
normalized  depth  function  at  source 
and  receiver  and  the  depth  function  and 
their  derivatives  at  each  layer  boundary) 
from  the  temporary  output  file  NFILE (2) . 
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If  NOPT  **  1 only  the  first  4 parameters 
are  needed.  NPPL (I ) gives  the  number 
of  data  points  in  each  layer  (the 
surface  is  not  counted  as  a data  point, 
data  points  at  layer  boundaries  belong 
to  the  preceeding  layer),  i.e.  if  XINCRl 
= 50  and  the  widths  of  the  first  and  sec- 
ond layer  are  100m  and  240m,  NPPL(l)  = 2 
and  NPPL (2)  = 5,  etc. 


Data  are  taken  at  intervals  of 
XINCRl  to  the  largest  multiple  of 
XINCRl  which  falls  into  the  layer  in 
which  XMAX1  lies.  ISLAY  is  the  total 
number  of  data  points  taken,  ISOU  and 
IRECi  specify  the  data  points  which 
corresponcto  the  source  and  receive r 
depths  respectively.  If  e.g.  the  source 
depth  is  100m,  the  receiver  depth  150m, 
in  the  example  given  above,  then 

ISLAY  =7  (=  NPPL (1)  + NPPL (2)) 

ISOU  = 2 
IRECI  = 3 


If  NOPT  = 2 the  depth  functions  a t the 
receiver  depth  are  obtained  in  a second 
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loop  through  the  programs  as  explained 
above,  if  the  receiver  position  does 
not  fall  in  an  up-layer  for  any  of  the 
modes  one  sets 

IRLAY  = IRLAX1  = total  number  of 
data  points  taken. 

If  the  receiver  depth  does  fall  into 
an  up-layer  for  some  modes  then  it 
is  necessarily  very  large  and  one  sets 
XMAX2  = total  vater  depth  H. 

Then 

IRLAY1  = total  number  of  data 
points  taken  (;£  H/XINCR2) 

and 


IRLAY  = IRLAY1  if  the  ocean  bottom 
is  not  a data  point, 

IRLAY  = IRLAY1— 1 if  the  ocean  bottom 
is  a data  point. 

(2)  Preparation  of  profile  data  file.  The  data  essentially  consist 
of  the  array  (1.2)  plus  the  corresponding  depth  values  Z^. 
Sometimes  the  sound  velocities  are  not  given  as  the  same  depth 
for  all  profiles.  Therefore  the  depth  values  are  also  taken 
as  a function  of  range. 
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The  profile  data  file  is  read  by  the  subroutine  SELAUT.  At 
present  the  following  data  structure  is  assumed: 

1st  line:  N = Number  of  layer  boundaries  for  first  profile^ 

FORMAT  (13)  . 

2nd-jth  lines:  N + 1 velocity  values,  5 to  a line,  FORMAT  (5D14.7) 
= sound  velocity  at  each  layer  boundary  + bottom  velocity 
for  first  profile. 

(j  + Dst-^th  line:  N + 1 depth  values,  5 to  a line,  FORMAT 

(5D14.7)  = depth  values  of  layer  boundaries  + bottom  depth 
for  first  profile. 

lst-kth  line  repeated  for  each  additional  profile. 

If  this  arrangement  is  found  inconvenient  it  can  easily  be 
changed.  The  only  program  affected  by  this  change  would  be  the 
short  subroutine  SELAUT. 

The  program  in  its  general  form  does  not  accept  isovelocity 
layers.  From  (3.6)  one  can  see  that  for  those  cases  the  Airy 
function  arguments  become  oo  and  (3.1)  shows  that  the  Airy 
functions  are  not  defined  there.  Isovelocity  layers  do  rarely 
occur  in  empirical  profiles  and  if  they  do  they  can  easily  be 
avoided  by  replacing  the  vanishing  slopes  by  very  small  slopes 
(e.g.  .01  •'m.  /sec  per  100m)  . If  desired  the  program  can 
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easily  be  generalized  to  include  the  (essentially  trivial  case) 
of  isovelocity  layers  by  adding  2 subroutines  which  propagate 
the  down-  and  up-  solutions  for  these  cases  between  layer 
boundaries  and  which  are  called  by  the  main  program  in  place 
of  PROP,  LARGE  or  START  and  PR.  or  LG. 

Before  setting  up  the  profile  data  file  it  is  necessary  to  check 

the  widths  of  the  vertical  layers  and  to  make  sure  they  do  not 

exceed  a certain  maximum  value  which  is  given  by  the  requirement 

that  the  Airy  functions  should  not  grow  or  decay  within  one 

layer  by  more  than  the  maximum  number  the  computer  can  handle. 

The  program  in  its  present  form  was  written  for  the  PDP10 

in  which  the  exponents  are  limited  to  % + 35.  At  the  end  of 

each  up-layer  the  depth  functions  are  scaled  down  to  10-^. 

Hence  the  maximum  allowable  growth  within  one  up-layer  is 

10^0  (or  1060  for  the  bottom  layer  where  the  depth  function  is 

-25 

made  to  start  out  with^lO  ) . Layers  which  do  not  satisfy 
this  requirement  have  to  be  subdivided.  An  estimate  for 
the  maximum  allowable  layer  width  can  be  obtained  by  using 
the  formulas  derived  in  section  3.  From  there  it  follows  that 
the  Airy  function  arguments  within  the  lowest  layer  (which 
is  usually  the  one  causing  trouble)  is 

— -A 


78 


= (2W2  a.A,)173  (vz>cmin- 

The  ratio  of  the  Airy  function  values  at  2 locations  within 
that  layer  is 

82  ~ e Vi 

A’  ‘^e  * *£*  (A,/«JV‘  «.Au)Kf  (^)!/*  - (^)V‘J 

where  /h  is  the  sound  speed  gradient  and  h the  width  of 
n n 

the  layer.  For  the  example  chosen  in  section  3 (f  — 250  UiL. 

A /h  = .015  1/3  and  h = 3000  m) : 
n n n 


Ai(Kz)  2.-Z3 

At(?;)  ~ 10 


4U 


' t - . 


With  Z = 3 = 4000  m,  Z.  = 1000  m as  in  the  example  one 

2 B JL 

obtains 

ArurJ  . 

: — \ ~ io 

A 

as  before  (see  section  3) . Hence  in  this  case  the  3000  m wide 
bottom  layer  would  have  to  be  divided  up  into  about  6 separate 
layers  of  **  500  m each. 

The  example  presented  is  actually  something  like  an  extreme 
case  sinfis  250  He.  for  a deep  ocean  case  is  about  the  upper  limit 
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of  application  for  the  normal  mode  theory  and  furthermore  the 
layer  thicknesses  as  they  are  usually  given  rarely  exceed  1000  m. 
Hence  only  in  rare  cases  do  layers  have  to  be  divided  up.  The 
program  presented  in  this  report  was  written  for  the  PDP10 
in  which  the  exponents  are  limited  to  a/  + 35.  For  a machine 
with  a larger  exponent  range  the  restrictions  on  the  layer  widths 
can  of  course  be  relaxed  correspondingly. 

(3)  Adjustment  of  array  dimensions. 

The  maximum  number  of  layers  in  any  of  the  profiles  determines 
the  dimension  of  the  arrays 

ISC,  MSC,  S,  B,  X,  Y,  T,  TD,  C,  P,  ALPH,  XS , XT,  NPPL. 
and  the  second  dimension  of 
NSC. 

I The  dimensions  of  these  quantities  should  be  chosen  as  the 

number  of  layers  + 2. 

The  dimension  of  the  arrays 
JUPA,  IRECA, 
the  first  dimension  of 
NSC 

and  the  second  dimensions  of 
UO,  UF,  XKA 

are  determined  by  the  maximum  number  of  modes  obtained  for  any 
particular  profile. 
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The  array 
XKST 

has  to  be  dimensioned  according  to  the  number  of  profiles. 
Dimensioned  arrays  are  limited  to  the  main  program,  except  for 
XKST  and  NFILE 

which  appear  both  in  the  main  program  and  in  the  subroutines 
PARAM  and 

X and  Y 

which  are  used  in  SELAUT. 

(4)  Execution  of  program. 

After  the  preparations  described  in  ( 1 ) — ( 3 ) are  taken  care  of 
the  program  can  be  executed.  The  following  routines  and  files 
are  involved: 

Main  program:  SEARCH 

Subroutines:  PARAM,  START,  PROP,  LARGE,  PR,  LG,  DAIRY, 

EXTREM,  SELAUT 

Profile  data  file:  NFILE  (1) 

The  results  of  the  calculations  is  output  through  the  files 
NFILE  (3):  depth  functions  at  source  depth 
NFILE  (4) : depth  functions  at  receiver  depth 
NFILE  (5) : eigenvalues 

NFILE  (6):  depth  functions  and  derivatives  at  layer  boundaries 
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(5) 


Mode  number  test. 

If  the  eigenvalue  spectrum  contains  almost  degenerate  modes 
it  is  possible  that  the  search  procedure  contained  in  the 
program  skipped  over  some  modes.  Therefore  after  execution  of 
the  program  one  has  to  carry  out  a check  to  see  whether  the 
eigenvalue  spectrum  for  each  of  the  profiles  is  complete. 

This  consists  of  rerunning  the  program  for  the  highest  mode, 
for  each  profile  (using  the  appropriate  starting  values  XKST 
for  the  eigenvalue  search)  with  a sufficiently  small  XINCR1 
to  obtain  a complete  picture  of  the  modes  and  plotting  the  modes 
as  a function  of  depth.  The  depth  function  values  can  be  taken 
from  the  temporary  file  NFILE (2) . If  one  finds  that  for  a 
particular  profile  the  number  of  modes  corresponds  to  the  mode 
number  then  one  can  conclude  that  the  solution  of  the  eigenvalue 
problem  was  successfully  completed  for  this  profile.  If  the 
number  of  modes  is  larger  than  the  mode  mumber  this  means  that 
some  modes  are  missing.  The  missing  modes  can  be  found  by 
repeating  the  procedure  described  above  for  some  lower  modes 
(e.g.  starting  with  a mode  number^  1/4  of  the  highest  mode 
number) , until  the  mode  number  of  the  missing  modes  has  been 
established.  The  exact  location  of  the  missing  modes  can 
then  be  found  by  running  the  program  between  the  neighboring  two 
modes  with  a sufficiently  small  XKD. 
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DIMENSION  S(60)|B(60).X(60)»Y(60)»T(60)|TD(60)«C(60) 

OIMENS  iON_Ei.60.li.  ALPHi-601».UO  ( 6...80.L.ilEi  6. 80-Lt.XkA(  6»U0 ) 

DIMENSION  XS(60),XT(60),NPPL(60),NSC(80,25),IRECA(a0) 

D4MENSlON_XKSia0UAirJLEi-ZJ 

COMMON  /IFLAGS/  MFLAG.MFLAGl.HTYPE  /JC/JCOUNT 

CQMMON_/SQNORM/ XNQRM 

COMMON  /DATCOM/XKST,  FRE0.RH01.RH02,xKD,XMAX1iXMAX2. 

_1 X.I  NCfil.AlJlCR2pAlPEUJM.LN0 . NJLLLEUXS  LAY  1 1RLA1 . 

2 ISOU. IREC.NOPT, IPRMIN. IPRMAX. ICT. IRLAYl, IREC1 

CALL— EAR  AM lt\ 

DM*U.0O0  ' r 

XKOS=XKO 

MTYPE*NFlLE(2> 

NtYPEsNFlLE(2J 

DO  15  IPR*IPRMIN. IPRMAX 

CLOSE — t UN  Ll=  NE-LL£-(-lJJ 

CLOSE  (UNIT cNF ILE ( 2 ) ) 

C ALUS E L A UJi lPJUNFlLE.a)  .-NFJLE(2 ) ) 

CLOSE  (UNITsNFILE(2)) 

XNQR1E0-..&Q0 

SCALE=1,0O0 

SCALE2;1i0D0 U) 

OSCAlEpI >0020 

LCOUN.T»0 

READ ( NF I LE ( 2 ) > 26 ) NUMBER 

F.QRMAL(J3J 

REA0(NF1LE(2).1)  (X(I),Y(I),I»1,NUMBER) 

REA0LNF_lLLL21.J.l_CaaiM 

CL0SE(UNIT=NFILE(2))  ' 

CMJN  = Y_a) 

00  199  I«l, NUMBER 


CONTINUE 

NMlN=NUMBER-j 

H»X(NUMBER> 

JD=NM I N"  JM 1 ND 

H2=H*H 

PI =3 , 1415926 53589783238400 
eoifVi,0D0 

_F  0_RM  Ali  20J 

0M=2 , 0O0“FREQ*P l 

OM  2*0  Mj»*2 

XK=0M/CMIN 

IF(XK.GT.XKSTUPR))  XK*XK$1 
P(l)i=H2#0M2/(Y(l)**2) 

00  43  l«i.NMIN 

X< I*1)=XC I *1 ) /H 
P ( I ♦! ) °H2 *0M2/ ( Y ( 1 *1 ) *»2 ) 
DO  44  iPl.NMIN 

-sxij=iPiitiupj.i  iimuixb 

C(!>=P(1>"S(I)*X(I> 

_CONTJ  NUE 

ITERa0 


I 


00560 

00570 

00500  

00590  66 

j)06a0_ 

00610 

00620 

00630  67 

00640  S90- 

00650 

00660 

00670 

00680 

00690 

00700  __ 

00710 

.00720 

00730 

00740 

00750  990 

00760 

00770 

0078.0 

00790 
00000  _ 

00310 

00820  _ 991_ 
00830 
00840  __ 
003^0  20”" 

00860  994 

00870 


00890 
00900  _ 
00910 
00920 
00930 
00940 
00950 
00960  _ 
00970 
00980 
00990 
01000 
01010  ' 
01020 
01030 
01040 
01050 

01060  

01070 
01080 
01090 
01100  _ 
01X10 
01X20  _ 
01130 

01140  

01150” 
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THs 1*00 0/3*000 1 

00  66  I >1 1 NM IN  I 

SN  = S ( t ) /DABS! SiXU , . 

ALPH(l)sSN»(S(I)*»TH)  V(h-) 

III!  =S  tlUOAESCStU) _ 

00  67  I s2 1 NN  I N 

Til)  s AlPH  U-ll/ALPHUJ 

CONTINUE  J f. 

IF  LAG-.0 — ...  71 

ICAL=0 

XHAX?XmAX1ZH 

X I NCR*X I NCR1/H 

iF  LLCQUHI.  *EILaCll_&a_Ifl_9L80 

MFLAG=0 

MF.LAGlrB 

XNORM=0,0O0 

JJER=0 

KFLAC50 

-XFACm*0O0 - 

XKsXK-XKQ  -ft 

U8QTMSRHQ2/BK01 ' 

IF(XK,GT. (OM/CBOTM))  GO  TO  991 

iE_LXKD  t-NEj.XKQ.Sl  GO  JO-  9 Eg 

XKsXK^XKD 

XKOSXKOS710 

GO  TO  990 

CONTINUE 

00  20  I=l»NMiN 

J SC.U  ).?.0 

MSC  < I ) =0  r 

. -CONTINUE ^ \ 

UK=OSQRT( (XK**2)-< (0M/CB0TM)*«2) J 

I F < NFL  AG . EQt-U-ANQ  J C 4U.N  E.dJ_J.U  EAjL  [ COUNT)=J  U P 

JTER=JTER*1 

X N = H 2JLIXK*  • 2.1 

00  3 J»1.NM1N  \ 

JUP  =0 : 

TTsT(J) 

IFJ  X ( J ) , G T , X H A.X.l_M FLA  G 1*0 

iF(J.EO.l)  XNEW  = -X I NCR 

SN  = S(  J)/PABS-(S(  JLL) 

AL=SN»(S< J)«*TH) 

Bt  J)SC(  JV»XN. 

BE»B<J)/S{j) 

DX*-Al»tX(J)»BE) 

0X1=-AL*(X( J*1)*BE> 

SWr_l..P0 - 

IF(J.Eo.l)  CALL  START(AL.BE,XU)#X<J*1>,U0UT.UP0UTi 

—1 XNEW  t x i NC  R , SW . TJ) i 

IFU.Eq.I)  CO  TO  812 

CONTINUE 

IF(J.LE.JO)  GO  TO  605 

I F1K  F L AG  ».E  Q_l3J___.G0 JQ—6.00 

I F ( DXl . LE , 7 ) GO  TO  605 

KFLAG»i 

J0*J-1 

GO  TO  600 

CONTINUE 

IF  ( OX  .GT  , 0M1_CALL_L  ARGEIALj-BE  . X ( J) , X ( J«1 ) .U1N.  _ 

1 U0UT,UP1N,UP0UT,XNEW,XINCR,TT) 
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01140  }r<Dx.ir.nM>  cap  propiai  .re.xi j) ,*t j*d fyiN,  ' 

01170  1 UOUT.UPIN.UPOUT.XNEW.XINCR.TT) 

01180  812  UIN=U0UT 

01190  UPlNsUPOUT 

01200  IF  (MFLAG.NF.l.OR.  ICAL.E0.1J  GO  TO  113 

01210  IFU.NE.I)  GO  TO  120 

0l2?0  CAl L OAIRYUAI  *BF.  AI  .AIP.BI  ,RJP> 

0 1 2 3 0 UOUT050, 

01240  UPOUT0=OABS(BI»AIP-AI*BIP) 

01250  JPL=1 

0l2ft0  WR I TE ( NF ILE  (_7J_»  130J  . JPL.KFLAG  ._X  ( L)_.  AL.  BEUIQUT0UJPOUT0  - 

01270  120  CONTINUE 

01280  JPLs J>1 

01290  WRITE(NFILE(7).130)  JPL , KFL AG » X ( J*1 > . AL . BE . UOUT , UPOUT 

01300  130  FORMAT ( 2 I 4 » 50  ) 

01310  115  CONTINUE 

01320  3 CONTINUE 

01330  SCALE=UBOTM/U!N 

01340  srAiF2c<ScALF*»2> 

01350  UPBOTM=»H«UK/'AL 

01360  FIINC  = U!N»I1PR0TM-UPIN*UB0TM 

01370  XNORM2=0 i 000 

01380  GO  TO  45  _____  _ _ . . . _____  ___ 

01390  600  CONTINUE 

01400  KFi  AG=1 

01*10  isw=icount-i 

01420  JX=J 

01*30  00  4 JsNMJN.JX.-l 

01*40  JUP=NMlN-JX*l 

31* ”0  FCT=1 i 000 

01*60  MS  = 0 _____  _____  _ ____  _ 

01*70  IF(MFLAGl.EQ.l)  MSsl 

01*80  I F ( X < J)  .GT.  XMAX.AN0.MS.E0.il  MFLAGH0 

01*90  I F < J.EO.NNIN)  XNEWsl.000 

01500  SN  = S( J)/OABS(S( J)  ) 

01510  AL=SN«(S< J)*»TH) 

01520  TOt J)SALPH( J*1)/ALPH( J) 

01530  IFU.tT.NHIN)  GO  TO  606 

01540  TO ( J ) =1 . 000/AL 

01550  UOINsUBOTm«XFACT/1,0O25 

01560  UPOlNs'-H°UK/-AL 

01570  UPDIN=UPOIN»XFACT/1.0O25 

01580  606  CONTINUE 



01590  B(J)sC(J)-XN 

01600  BE  = BU)/S(J7 

01*10  DX1=-AL*<X( J)*BE> 

01620  0X=-AL«(X(J*1)*BEI 

01*30  jF(OX.GT.DM)  CALL  LG ( At ,BE , X U*1 ) , X U) . UDJ N, 

01*40  1 UOOUT  lUPOlN.UPOOUT  . XNEW.XJNCR.  TQU)  »XS(J) . XT  ( J i ) 

01*50  IF(DX.LE.OM)  CALL  PRCAL.BE.X(J*1>iX(J).U0IN, 

01*60  1 UOOUT.  UPOIN,UPOOUT.XNEW.XINCR.TOU>.XS<J),XTU 

lit 

01*70  IFU.EO.JX)  GO  TO  300 

01*80  IF(MFLAG.E0,1 ) GO  TO  110 

01*90  JF< 1FLAG.E0.1)  GO  TO  608 

01700  603  !F< IFLAG.NE.l, AND. OABStUDOUT) .GT. 1.0-3) I SC  < J ) 

I sc  < j > *: 

01^10  I F ( I SC ( J) ,EO.0)  GO  TO  300 

01720  UDOUTSUOOUT/I.010 

01730  UPOOUTsUPOOUT/1,010 

01740  IF (OABS(UOOUT) .GT.l.D-5)  GO  TO  603 

01750  GO  TO  300  J 
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01760  _ .608 
01770 

01780 

01790  110 

..Jl 1.000 

01010  604 

01020 

01030 

.01040. 

01050 

01060 

01070  301 

01080 

01090  300 

01900 

01910 

_ 01?20 

01930 

.01940 

01950 

.01960 

01970 

...01?M 

01«90 
02000  4 

02010 
02020 
02030 

02040.  . 

02050 

02060 

02070 

02080  _28  _ 
02090 
02100 
02110 
02120 

02130  46 

02140  _ 47 
02150  22 

02160 
02170 
02180_. 

02190 

02200 __  _45_ 
02210 

02220 _ C_ 
02230 

02240  _ 

02250 

02260  

02270 

02280 

02290  "C 

02300  _ 

02310  170 

_02320 

02330 

02340 

02350 


.-UOOUT*UDOUT/10..D0«*C-ISCUU10J A ' 

UPOOUT=UPOOUT/10.O0««< ISC< J)*10) 

GQ  T 0.300 

CONTINUE 

IfLtIC  AL.  EQ  tl) — G Q_I.O—301 

IF  (OABS (UDOUT ) ,GT, 1,0-5 )MSC<J)=MSC<J)«l  > (lo\ 

I F( MSC  CJ ) , EQ .0  J_ GO  -TQ.-300 

UDOUT=UDOUT/1.O10 

UPOQUTsUPOQUX/1,0.10 

IF(OABSIUOOUT) ,GT, 1,0-5)  GO  TO  604 

GQ_T.O_ 300 

UDOUT  = UDOUT/10.D0»*(MSC(J)*10)  M'i) 

UPDOUT;  UPOOUJ/10 . 00 * IMS  C.UI«10 1 

CONTINUE  J 

JP.LRJil 

UX  = 0 

PX?0, 

1FU.EQ,  JX.AND.  JX.NF.NMIN)  UXsUDJN 

iF-(d . EQ . JX.  AND . JX»NE_.NM  IN1_PX*UE0IN 

IF(ICAL.NE.l.ANO.MFLAG.EQ.l)WRITE(NFiLE(7),130> 

1 JPLiKFUG.  XUJj-AU5E..UX.P-X 

UOIN=UOOUT 

UPDlNrUP.OOUI 

IF (MStEOtl)  MFLAG13! 

CONTINUE - 

UDOUTsUD OUT/1, 0015 

_ . UPOOUTeUPOOUT/l  ,0015. 

FUNCsUlN«UPOOUT-UPlN*UOOUT 

IF  LMFUG . NE-.-11.GQ.  _T0_45 

NS  = 0 

IF(JX,EO.NMIN)  G0-TQ_22 

00  28  I 3 JX*1 , NM IN 

NS  = NS-MSC(I1*ISC.(JD 

00  47  K*JX.NMIN  , . 

-._nsc<icount.ki=ns yjlU 

IF(K.LE.JX)  GO  TO  47 

DO  46  NeJX,K 

nsc  (i  count, k)=nsc< icount,k>*msccn) 

CONTINUE 

XNORM2=0,D0 

XST*XS(  JX)*XT  l JX). 

IFINS.EO.0)  XN0RM2sXST 

IF « NS , NE , 0 , AND,  DLOG10  ( XSD.-20 ,.O0*NStGJLuHL30.,O01 / 

1 XNORm2=10,O0««(OLOG10(XST)-NS*20,O0>  J 

CONTINUE 

!F< ICAl.EQ.I)  GO  TO  892 

!.YPE_1,XK,F.UNC 

XNOR1S0,0D0 

IF (KFUaG,EQ.0)  XNORl  = RHO2/{2..0b0»lL«l!kJ 

XNR  = XN0R1*(XN0RM2*XN0RM)«SCALE2  (TTi-v 

Iff XN.R , G T.0.000)  SON3<1.0O0/DSQRT(XNR))*DABS(SCALE)  AM* 

IF (MFLAG.EQ.l)  WRlTE(NTYPE.l)  SQN  J 

_ ..I  F.<  NFt AO , NE , 1 )__G0— TQ-_1_65 

STORE  U.Up  AT  LAYER  BOUNOARJES  FOR  MODE  JCOUNT 

WR  ME  f 21 , 170J— NU.NBER 

FORMAT< 14) 

CLOSE  I UN  I T£29J 

00  160  INsl.NUMBER 

RE AOJNFJ L E ( 7), 130)  JC ,KC.XC,ALCtBEC,U*C,UPXC 

UXC»UXc»SON 


I 1 1 .'l 
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02360  UPXCfUPXCaSQN 


C 02370  I F <KC • EQ • 1 ) UXC»UXC»10*»(-10*N$C( ICOUNT, JX> ) 

02380  IF(KC.EQil)  UPXC-UPXC«10**  < -10*NSC { I COUNT  1 JX  > ) 

_ 02390  160  HR  I TE ( 21 1 132 ) XC . ALC . BEC . UXC . UPXC 

C 02400  132  FORMAT (50) 

02<10  165  CONTINUE 

02»20  CLOSE ( UN  I T = NF_ILE_(_7 ))..._  . 

C 02430  IF(MFLAG.NE.1.0R.N0PT.E0.1>  GO  TO  892 

02440  I CALC1  ...  _ _ 

02450  XMAX=XMAX2/H 

C 02460  XINCRSXINCR2/H 

- (13) 

02470  mflagi=i 

02480  GO  TO  994  _ . 

- /,  \_ 

< 02490  892  IF(MFLAG.EO.l)  GO  To  890  < {IV 

02500  IF(OABS(EDIF)1LTJ.1.0D-15)  ITTR»1 1 

02510  !F(ITER,EQ,0)  GO  TO  444 

< 02520  ICnUNTslCOUNT+l 

02530  EOlF=1.0D0 

02540  MFLAGlnl 

f 02550  MFLAG*1 

02560  . IF_<KFl,AC.EOa)  X T. ACJZU.IN/.CiiQ QUIP  1 13 Dl5J 

l('s) 

, 02570  WRITE(NTYPE,1)  XK 

t 02580  TYPE  770i XK^LCOUNT 

02590  tTER=0 

02600  GO  TO  994 

C 02610  444  CONTINUE 

02620  FU(JTER)=FUNC  " 

02630  EIGIJTERJsXK 

C 07640  IF«IFL*G.F0.1  ) fifl  TO  880 

- 026^0  I F < JTER  , EO , 1 ) GO  TO  990 

- 02660  PRO0  = Fu<JTER)/DABS<FU( JTER)  ) 

C 02670  PR0D=PR00»FUC JTER-1)/DABS<FU( JTER-l) ) 

02680  IF  J PROD  t GT .0.0OB)  GO  TO  990 

. 02690  iflag=i 

* _ 02700  XK  = E I C { JTER.)  -DSQRT  ( F_UC  JTLR ) /_l£m  JIEB  J "Fll.tjTER.-Jl 

... 

02710  1 *(EIGIJTER)-EIG(JTER-1)  ) 

02720  GO  TO  994 

(it) 

C 02730  880  JMINf JTER-l 

02740  KCJMIN 

rvbl 

02750  GO  TO  882 

C 02760  00  881  K=JMIN.i.-l 

02770  PR0D=Fu< JTER) /DABS (FU( JTER)} 

_ 02780  _ . PROD=PROD*FUIK)/OABSirULKLL_ 

C 02790  IF(PROo.LT.0.0D0>  GO  TO  882 

02800  881  CONTINUE 

02810  882  SE  = {FUUTER)-FU<K))/(EIG(JTER)-EJC(K)) 

C 02820  SIN*FU( JTER)-(SE»EIG( JTER)) 

\ 

02830  XK=-S JN/SE 

02840  EOIFbXk-EIGIJTERJ 

; 

1 

C 02850  770  FORMAT (iO» 15) 

02860  GO  TO  994 

1 

ti  02870  980  CONTINUE 

C »i 02880 J_!P E _3 1 # J J U p A (J )_». !_■  Li  I C 0 U Nil 


— — 1 j (U.J  

to  02890  31  FORMAT ( • UP-LAYERS  **/10I4> 


9 

02^00 

CLOSE  (UNIT«NTYPE) 

- - 1 4 

c 

6 

02910 

NTTF0 

7 

02920 

DO  41  I>l,NMiN 

6 

02930  41 

NTTbNTT*  NPPL ( I ) 

c 

* 

02940 

00  16  I »1 f I COUNT 

4 

02950 

IF < I . GE « I SW , OR , NOPT , EO . 1 ) IRL*YbIRLAY1 

J 

3 


c 


( 

( 

t 

( 

t 

( 

( 

(. 
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{!?960  f RpC  = NTT"  IRf!''l  + l 

02970  DO  42  jsl ,'  NM  I N- JUPA  < I ) 

02980  42  IREC3 lREC*NPPL(J) 

02990  95  I P ( I . GE , I SW)  IREC=IREC1 

03000  1 RECA I l ) 5 1REC 

03010  READ(NTYPE.l)  XKA(IPR.I) 

03020  no  70  Ksl(  !<?|  ay 

03030  READ ( NT  YPE • 2 ) BXTMP 

03040  i F t K . En . I sou > RXsRXTMP 

03050  70  IRONED, IREC. ANO.NOPT.EO.l)  CX»BXTMP 

03040  RrAD(NTYPr.t>  YN 

03070  lF(NOPT.EQ.l)  GO  TO  30 

03080  DO  71  Ksl, IRLAY  . 

03090  READ ( NT  YPE  « 2 ) CXTMP 

03100  71  IF  (K.EO.  IREC)  CX=CXTMP 

03110  30  CONTINUE 

03120  NSPPLs0 

03130  DO  39  Lsl.NMIN 

03140  NSPPL^NSPPL  + NPPLILl - - - . . .. 

ST)  _ \ 

03150  Kl=L 

03160  IF (NSPPLiGT . IREC1)  GO  TO  40 

03170  39  CONTINUE 

03180  40  CONTINUE 

03190  UO(  IPR.  I >=F3X»XN 

03200  16  UF(IPR.n=CX*XN*10,D0"*E“10*NSC(liKLU  ..  - 

03210  CLOSE  (UNJTeNTYPE) 

03220  19  XKO*XKDS  - 

03230  TYPE  4e, ( IRECAI I ) , 1=1. ICOUNT) 

03240  48  FQRMATC  u OF  REC  IN  WFCT  FILE  3 '/10I4) 

03250  TYPE  18,  IPR.  ICO’.'NT,  ISW 

03260  18  FORMATC  END  OF  PROFILE  J* ',  .1.2 _ * OF  EV  I 

N THIS  P 

03270  IROFILE  =',!2.'  I SWS ’ ,12) 

03280  15  CONTINUE 

03290  DO  17  jslPRMIN, IPRMAX 

03300  WR I TE  < NF  I LE ( 3 ) , 7 ) ( UO ( J , I ) , I =1 , I M ) 

03310  WR I TE  « NF I LE  C 4 ) ,7)  < UF ( J , I ) , I =1 , J M > 

03320  WRITE<NFILE<5),7J_  (XKAlJ,JUl»liIJ4J.  ..  . 

03330  17  CONTINUE 

03340  2 FORMAT  ! 020 . 9 ) 

03350  7 FORMAT ( 860 ) J 

03360  996  STOP 

03370  END 

1 
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00010 

00030 

00030 

-00040- 


00050 

00060-  _10 
00070 

,00080 

00090  20 

-00100 


00110  1 
-00120 Z- 


SU9R0UTINE  SEUAUT<  JP.NFIL.NTEHP) 

DOUBLE  PrtECI  S I ON  -X  f55  X r-Y  ( 55  X 

00  10  I«?1.1P 

-BE-Aaoip-twax-N 


READ(NriL.l)  <Y( J) »J«1.N*1> 

- HE  A0(  NP  1 L 1 1 >-  ( X { >>  r J=  l^N«-tX- 
WR I TE  < NTEMP ,2)  N 
-00— 20-  I si  ^N- 


WRITE(NTEMP,3> 
-WRUttNTCHPrl-X- 
FORNAT (5014,7) 
-£OmU*-(-l-3->- 


X( I) ,Y< I) 


00130 

-00140- 


FORNAT ( 2D ) 
-RE-T-ORN 


00150 


END 


( 

( 

( 

< 

( 

< 


c 

( 

< 

< 
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00010 
00020- 
00030 
-20040- 
00050 
. 00060  - 
00070 
-000(10- 
00090 
—00100— 
00110 
-00120- 
00130 
-00140- 


SUBROUtINE  param 
-IMPLICIT  DOUBLE  PRECISlOM-tA-H^)-24- 


0 IMENS I ON  NPPL<60) »XKST<10) ,NFILE<7) 

-COMMON-/ D AlcflH/  XK SX^PTro^RMQ-l-^RMOZ > X (CO t XMAXl-^XHA 42  i 


1 XINCR1,XINCR2.NPPL. JMIMD.NFILE. ISLAYi IRLAY* 

-2  -I  SOU  r I BEC  , NOPT - I PRKIW.4PRHAX.-I  C-T-»-IBL  AY-1  ,-IRECl- 
OATA  NPPL/10,14,24,57*0/ 

OAT  A-X  K ST  / 1 0 a-l  000  .-00/- 


00150 

-00160. 


DATA  MFILE/38,24.25,27.28,21.29/ 

4^=^043.7804 

FREQ=5.D0 

_r KOI- 1^00 

RH02=1.00 

_XHAJO-=-l  00-. 

XINCR1=50, 

-XMAX2  = 10g, 


00170 
-001 00- 
00190 
-00200- 
00210 
-00020- 


JSLAY=10 

-IOLAY55-1- 


IS0U=1 
-1RECU2- 


IRLAY1=51 


00230 

-0024.0- 

00250 

-00260- 

00270 

-00260- 

00290 

-00300-- 

003lf 


XINCR2=1O0. 
XKO=-4d-4— 
ICT*20 
JHINO-10 


I PRM I Ns  1 
-jRRHA-Xs-A- 
NOPT  = 1 

RETURN 

END 
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00010 

00020 

00030 

00040 

00050 

00060 

00070 

00080 

00090 

00100 

00110 

00120 

00130 

00140 

00150 

00160 

00170 

00180 

00190 

00200 

00210 

00220 

00230 

00240 

00250 

00260 

00270 

00280 

00290 

20-00 

003i0 

00320 

00330 

00340 

00350 

00360 

00370 

00380 

00390 


C FOR 
C 

C 

c 


00410 

00420 

00430 

00440 

00450 

00460 

00470 

00480 

00490 

00500 

00510 

00520 

00530 

00540 

00550 

00560 

00570 

00580 

00590 

00600 


subroutine  oairyiox. ai.aip.bi.bip) 

IMPLICIT  DOUBLE  PRECISION  (A-H.0-2) 

DOUBLE  PRECISION_ARf.UMENTSt_T.HlS  ROUT i.NE_C ALCULAIES._THE._A 
function  aj<x>  and  its  derivative  aipui,  it  also  finds 

The  other  REAL_LINEARU  INDEPENDENT.  S0lU110N_BUX)_AND 

its  derivative  bipix) . 

the  definitions  ano  normalizations  are  as  in  .nbs  handbook- 

OF  MATHEMATICAL  FUNCTIONS, P, 446 

the  methods  used  are  power  series  .expans ion_f.or_small_x 

And  GAUSSIAN  integration  FOR  LARGE  x 

DIMENSION  X ( 16 ) , W ( 16 ) , XSQ (16 ) 

double  precision  dx,ai.aip,bi,bjp 

double  precision xs  ,xcube,aisum,aipsum_ 

oouble  precision  of.ofp.og.ogp 

Double  precision .f jm2 . EjMUF.j.tr jpju f.jp.2, factor 

double  precision  ci,c2. roots 

double  precision  dzeja.oarg.drootx 

double  precision  root4x,s,co. ratio. efac.hetaso 

double  precision  sumr,sumi.sumrp,sumip,termr,termj 

double  precision  dzero.oa.ob.oen.one 

Double  precision  x,w,xsq 

double  precision  rsq.temp,  rtpi,rtpi2 

double  precision,  terma, terms 

logical  needbi 

oata  dzero,one  /0, 000, 1,000/.. 

Data  R00T3/1. 73205080756887700/ 

DATA  Cl  * C2  / • 35502 8053 887 81 7D0,._»2508194037.92807D0/_ 

DATA  RTP I /, 282094791773878100/ 

DATA  RTP  12/.  56418958354  7 756200/ 

itigns  and  weights  for  ib-term  sum  for  airy  functions 

DATA  W(  1)  / 3.1542515762964787D-14/ 

DATA  W<  2)  / 6.63942108195849210-11/ 

DATA  W<  3 > / 1 • 758388 906 134 56690" 08/_ 

OATA  W ( 4)  / 1.37123923704358150-06/ 

DATA.WI  5 ) J 4,43509666392843500-05/ 

OATA  W(  6)  / 7.15550109177182550-04/ 

DATA  W?  7)  / 6.48895661033353810-03/ 

DATA  W(  8)  / 3.6440415875773282D-02/ 

DATA  VL  9>_/_1, 43997924105909990-01/ 

OATA  W < 10  > / 8.12311413362614860-01/ 

OATA  X<  1>  /_JL,  40830010721809640*01/ 

OATA  X<  2)  / 1.02148854791973310*01/ 

DATA  XI  3)  / 7, 44160184504509300*00/ 

OATA  X(  4)  / 5.30709430617819270*00/ 

DATA  XI  5)  / 3.63401350291324620*00/ 

DATA  XI  6)  / 2.33106523030524500*00/ 

DATA  X(  7)  / 1,34479708246092680^00/ 

OATA  X<  8)  / 6,41888583695672960-01/ 

DATA. .XI  9)  /__2, 01003459981210460-01/ 

OATA  X<10)  / 8,05943591720528330-03/ 

0ATA_XSQ(_1 1_./ 0,198333172485621700-03/ 

OATA  XSQ ( 2)  /0, 104343885353116500  03/ 

OATA  XSQl_3j„/0, 553774380201781700  02/ 

DATA  XSQ ( 4)  /0.28165249974668990O  02/ 

DATA  XSQ < 5)  /0, 132060541393558000  02/ 

OATA  XSQ ( 6)  /0, 543386510793804400  01/ 

data  xsqi_7>  /0, 180847919299542000  01/ 

DATA  XSQ(  8)  /0. 412020953878836900  00/ 

DATA  XSQ ( 9)  /0 , 404023909244180700-01/ 

OATA  XSQ( 10 ) /0, 64954507 303538390D-0 4/ 


2623JL 


OATA  X<  II  / 
OATA  X<  2)  / 
DATA  XI  3)  / 
OATA  X<  4)  / 
DATA  XI  5)  / 
DATA  XI  6)  / 
DATA„X(  7) J 
DATA  XI  8)  / 
DATA  XI  9)  / 
OATA  X<10)  / 


26250 

26260U 

26270 

26280L 

26290 

.26300- 

26310 

.26320- 

26330 

26340- 

26350 

26360- 

26370 

26380- 

26390 

26400- 

26410 

26420_ 

26430 

26440L 

26450 

26460L 

26470 

26480L 

26490 

26500_ 

26510 

26520_ 

26530 

26540_ 

26550 

26560_ 

26570 

26580- 

26590 

26600- 

26610 

26620- 

26630 

26640_ 

26650 

26660_. 

26670 

26680 

26690 

26700- 

26710 

.26720- 

26730 

26740 

26750 

26760- 

26770 

_26?60- 

26790 

26800- 

26610 
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00610 

C. 

PQSI_T.LQNS  AnQ-UE-IERTS.  FOR 4-TERM  Sl,IM  Fl)R  4 I RY  FUNCTION* 

26820 

f 

00620 

DATA  W<11)  / A.77639030575772630-05/ 

26830 

r 

00630 

DATA  W(l2)  / .4.99143064329109590-03/  _ 

26840  ... 

00640 

DATA  W(13>  / 8.6169846993840312D-02/ 

26890 

I 

.00650 

OATA.  Wll4J  Z _9.*0S79095flA5981  1090-01/ 

26R60 

f 

00660 

DATA  X( ll ) / 3.91983295544550916*00/ 

26870 

00670- 

DATA  X ( 1 2 ) _ /_ 1.69156190048235040*007 

96080 

■ 

00680 

OATA  Xll3>  / 5.02755324672630180-01/ 

26090 

j 

00690 

. 0ATA  X(14)  Z 1_.  924  70605620 156920-02  7 

26900 

I 

00700 

DATA  XSQ(ll)  /0. 153650903985966700  02/ 

26910 

.00710- 

- Data  XSQ(l2l-/0. 2861381.66316346100  Ot/ 

26920 

00720 

OATA  XSQ(13)  /0. 252762916486681800  00/ 

26930 

I 

00730 

DATA_XSQ(l4)  Z0. 37044934027789980(1-03/ 

26940 

i 

00740 

C 

POSITIONS  And  weights  FOR  2-term  SUM  FOR  AIRY  FUNCTIONS 

26950 

i 

00750 

_ DATA  W(l5)._Z  . 9. 6807280595  7 736040-01/ 

26960 

L 

00760 

OATA  W(16I  / 3 • 1 92 7 19 40 4 639580-02/ 

26970 

I 

_007.70 

OATA  X(l5)  / 3.68006018661530440-02/ 

26980 

| 

00780 

OATA  X(16)  / 1.05924693821123780*00/ 

26990 

00790 

DATA  XSQ(l5)  /0d3542842.977111070O-02/ 

27000 

i 

00800 

OATA  XSQ < 16 ) /0, 112200407610988100  01/ 

27010 

j 

0080? 

IF (OX. GT, -1000. 001  GO  TO  991 

1 

00804 

CALL  SlNC0S(0X.AI.A!P,8l.BIP) 

I 

00806 

RETURN- 

i 

00808 

991  CONTINUE 

1 

00810 

IFIDX.LT, -5,000)  GO  TO  100.. 

27020 

I 

00820 

NEEDBI=»FALSE.  ” “ 

27030 

t 

00830 

IF (OX >GT . 3 * 700 ) GO  TO  200  . 

27040 

i 

00840 

C 

this  route  for  smallx,  using  power  series, 

27050 

s 

00850 

C 

initialize  _ - ... 

27060 

i 

00860 

10 

xs  * ox*ox 

27070 

< 

< 

00870 

XCUBE  a XS  »DX  ..... 

27080 

s 

00880 

XS  = XS  *0.500 

27090 

} 

O089O 

OF  a Cl 

27100 

L 

00900 

OFP  = Cl*XS 

27110 

3 

00910 

DG  a C2*0X 

27120 

j 

00920 

OGP  = C2 

27130 

1 

00930 

A I SUM  a OF  - OG 

27140 

9 

00940 

AIPSUM  = OFp  - OGP 

27150 

j 

00950 

Bl  = OF  ♦ OG 

27160 

i 

00960 

BIP  s OFp  ♦ OGP 

27170 

i 

00970 

FJM2=-2.0D0 

27180 

| 

00980 

20 

FjM2=FjM2*3,0O0 

27190 

) 

00990 

FjMl=FjM2*0NE  ... 

27200 

9 

01000 

FjaFJMl*ONE 

27210 

I 

01010 

FJP1=FJ*0NE 

27220 

1 

01020 

FjP2=FjPl*0NE 

27230 

01030 

ratio  = xcube/fj 

27240 

01040 

Of  a 0F*RATI0/FJM1 

27250 

01050 

OFP  a DFP«RaT IO/F JP2 

27260 

01060 

DC  a DG#RAT IO/F jpi 

27270 

01070 

DGP  » DGp*RATIO/FJMg 

27280 

12 

01080 

Bl  a Bl  * <DF*OG> 

27290 

11 

01090 

BIP  a BIP  ♦ <DFP*OCP)  _ 

27300 

10 

01100 

I F < NEEDS  t ) CO  TO  80 

27310 

l 

0 

OHIO 

A[SUM  « AJSUM  * (OF-OC) 

27320 

1 

e 

01120 

AIPSUM  = AIPSUM  ♦ (OFP-OGP) 

27330 

•y 

01130 

C 

convergence  test 

27340 

6 

01140 

80 

IFIOABS(DF) ,GT. 1.00-16)  GO  TO  20 

27350 

6 

01150 

C 

convergence,  compute  functions 

27360 

4 

3 

01160 

99 

Bl  a RO0T3#8 I 

27370 

3-3 
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0U70 _ BIP  = R00T3#BIP. 

01180  C THIS  RETURNS  IF  X IS  BETWEEN  3.7  AND  8.0.  SINCE  (N  SUCH  CASES  MORE 
0U90...  C.  ACCURATE  VAlU?S_of.  Al_JANQ_A.If!_HAV.t_ALBEADjL.B£EH..£OUNO_BJLCAUSSlAN_ 
01200  C INTEGRATION 

01?i0_ IF  tNEEDB I ) RE-tURN 

01220  A I > AISUM 

01230 ._  A IP  = AlPSUH 

01240  RETURN 

01250  C GAUSS  I AN  . INTEGRATION  .F.Ofl_.LARGE_NEGA±LVE_X . 

01260  100  OROOTX  = OSQRT(-OX)  * 

01070 R00T.4X  .*  ..OSqRUORQOTX) 

01280  OZETA  » -.6666666666666667«DX*DR00TX 

01090 DARG.s  0 ZEI A—- _ *7.8 539816339 24  4.8.3 

01300  SUMR  = DZERO 

01310 S U Hi.  _5  . Q.Z  E R 0 

01320  SUMRP  * OZERO 

.01330 SUM  IP..?  OZERQ 

01340  C TEST  TO  SEE  HOW  MANY  TERMS  ARE  NEEDED  IN  GAUSSIAN  INTEGRATION 

01350  . IF10X.LT. {'200.00))  .GO  TO-140 

01360  lFtDX.LT. (-15,00) ) GO  TO  130 

0137.0 C IHI5  CASE  ..EoR.OX.BEIWEEN_r5  .0_ANO_ri5.t0 ; 

01380  LIML0=1 

01390 LIMH1  = 10. 

01400  Go  TO  149 

01<10  C THIS  CASE.F0R._DX  BET_WEEN_-15 ,0_AND._?200., 

01420  130  L1ML0=H 

01430  LIMHI=14 l! 

01440  GO  TO  149 

01450  C THIS  CASE  FOR  OX.LT-o*20i. 

01*60  t40  L I ML03l5 

01470  LlMHIsl6 ! 

01480  149  ZETASQ=0ZETA««2 

01490  ....  DO  150  K = UMLO.LIMHI_ 

01500  TERMR=W(K)X( (ZETAS0*XS0(K))**2) 

01?10  . _.SUMR  = SUMR  ♦-IERMR 

01520  TERMR=TERMR»X<K) 

01530  Sum i =sum i ♦Termr 

01540  TeRMR=TERMR«X(K) 

01550  SumRP  = SUmP  PtlER  MR 

01560  150  SUMIP=SUMJP*TERMR»X(K) 

0 1 570 SUMR=  ( S UM  R * ZE  T A S Q ♦. S UMRP)»ZETA  SO 

01580  TEMP=SUMI#ZETASQ 

.01590 SuMI  = <HMPtSUMJR)J*.DZ£IA ! 

01600  SUMRP=SUMRP*OZETA 

01610  SUMIP*SUM1P»TEMP 

01620  C FORM  AIRY  FUNCTIONS 

01630  196 S _s  OSIN(QARG) ' 

01640  CO  = DCOSWRG) 

01650  RATIO  J*  RTP.I2/R00T4X 

01660  A I = RAT  10* ( CO*SUMR  * S»SUMj)  • 

01*70 Bl  .=_  R AT  1 q» ( CO»SUM  1 - S»SUMR) 

01680  SUMRP=SUMRP*SUMRP 

01*90 RAIIO.  » - . 25D0/DX 

01^00  factor  * -rtpi2*root4x 

01710 _AJP  s^RATIO*AJ_r_  DROOTX*BI_t_FACTOR*ICO«SUMRP*S*SUMIPJ 

01720  BlP  s RAT  1 0*B I ♦ pROOTX»AI  ♦ FACTOR*<CO*SUMtP-S*SUMR>) 

0i730_ return , 

01740  C GAUSSIAN  INTEGRATION  FOR  LARGE  POSITIVE  X 

01750 200  OROOTX  *_PSQRT(DX1 

01760  OZETA  ■ .6666666666666667*DX»0RQ0TX 
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27510 

27520 

27531 

27  550 
27560 

27571 
_ 27580 

27590 
27600 

27610 

27620 

27630 

27640 
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. .27660 

27670 
_ -27680 
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27710 

27720 
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La 

01780 

01790 

RQ0T4X  = OSQRTCOROOTX) 

A?  = 02ER0 

27990 

.28000, 

< 

01800 

01010 

8 I = OZErO 

A1P  DZERQ 

28010 

28020. 

01820 

01830 

B IP  s OZERO 

IP(OX.LT.8.0D0)  . NEEOBI*.-TRUE.  . 

28030 

28040 

< 

01840 

01050 

C 

TEST  TO  SEE  HOW  MANY  TERMS  ARE  NEEDED  tN  GAUSSIAN 
JF(DX»GT. 15.000)  GO  .TO—230 

INTEGRATION 

28050 

-28060. 

< 

01860 

01870 

C 

THIS  CASE  FOR  OX  BETWEEN  3.7  AND  15, 

LIML0*1  - 

28070 

28080. 

01080 

01090 

LIMHi»10 

GO  TO  249  _ . ..  . 

2B090 

28100. 

( 

01900 

01910 

C 

230 

THIS  CASE  FOR  OX  GREATER  THAN  15, 

L I ML0*1 1 

28110 
. 28120. 

\ 

01920 

ei930 

249 

LIMHI»14 

00  250  KsLIMLO.LIMIU- 

28130 

28140 

01940 

01950 

0A=0ZETA*X(K) 

terma  « WlK)/QA.  . 

28150 

28160 

01960 

01970 

AI  * AI  ♦ TERMA 

A (PsAlP*TERMA*X(K)/OA 

28170 

28180 

V 

01980 

01990 

I F ( NEEDB 1 ) GO  TO  250 

0R=0?ETA-X(K) 

28190 

28200 

02000 

02010 

terms  s w(K)/ob 

Bj  = si  ♦ terms  _ __  _ .. 

28210 

23220 

{ 

02028 

02030 

250 

BIP=BIP*TEHM8*X(K)/DB 

continue 

28230 

28240 

02040 

02050 

C 

FORM  FUNCTIONS 

FACTORaRTP  J »0ZETA/R00T4X_.  ..  _.  . 

28250 

28260 

02  ?60 
02070 

Ratio  = 0,25Da/ox 

aIsAI»EFaC*FACTOR  - .... 

28270 

.28280. 

* 

02080 

02090 

C 

AIP=-(OROOTx*RATIO)»AI*RTPI«ROOT4X*EFAC*AIP 

THIS  IS  SATISFIED  ONLY  FOR  X BETWEEN  3 . 7 AND  8.0.  _ 

IN  THESE  CASES _ 

28290 

28300. 

02100 

02H0 

C 

C 

THE  01  AN0  RIP  ABOUT  TO  BE  COMPUTEO  ARE  NOT  SUFFICIENTLY  ACCURATE, 

THUS  RETURN  TO  POWER  SERIES  FOR  BI  AND  BIP. 

28310 

28320 

02120 

02130 

IF(NEEDBI)  go  to  10 

FACTORsFaCTOR+FACTOR 

28330 

28340 

( 

02140 

02150 

B1=9I*FAcT0R/EFAC 

B IP  = (OR0OTX-.RAT  10)  *B1-RTP  lE^ROQT-AXlB  IP/EFAC 

28350 

28360 

( 

32160 

02170 

RETURN 

end 

28370 

28380 

( 

( 
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00010 

00020 

00030 

SUBROUT  INF  E«TBFM(*r*f*PfB*BPi?FTA) 

1 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
niMENSJON  X 1 ( 4 ) i U ( 4 ) 

1 

00032 

00034 

00036 

00030 

00039 

00040 

IF  ( X iGT  1 15 , 00  > GO  TO  100 

CALL  OA!«Y<X.ArAP,n,nP' 

ZETA>1.O0 

RETURN 

■ 

( 

100 

CONTINUE 

VI  t lias,  9lQRS?95S44«5509in0 

( 

00050 

00060 

XI <2>*i, 691561900482350400 

XI  (3)=5,0?7553?4  672630ianr.01 

00070 

00080 

X I « 4) =1,92470605620156920-02 

U< 1 1*4.77639030575772630-05 

( 

00090 

00100 

W(2)»4 .99143064329109590-03 

U<3>=8. 61698469938403120-02 

( 

00110 

00120 

WI4) *9,08790958459811020-01 
nR00TXsnSQRT(x> 

00130 

00140 

ZETA=0.666666666666667O0*X*OROOTX 

RnnT4X=oSoRT^nRO0T^> 

( 

00150 

00160 

AS0.0O0 

APS0.0D0 

( 

00170 
00180. . . 

8*0.000 

RP=0.0n0 

00190 

00200 

00  1 K-1.4 

0As2ET**Xl Ck» 

( 

00210 

00220 

TERMA=W(K)/OA 

AsA*TERMA 

f 

00230 

00240 

AP=AP»TERmA«XI (K)/0A 

OB=2ET  A»X f 1 

14  1 l.i 

00250 

00260 

TERM8=W<K)/0B 

B*8*TERMB  . . . . 

K 

00270 

00280 

1 

BP=BP+TERMB*XI <K)/DB 

CONTINUE 

i 

00290 

00300 

RTP I =0.282094791773878100 

R TP  1 2 ’0.56418958354 7756200 

00310 

00320 

FAC=RTpi*2ETA/R00T4X 

RAT=0 • 25000/X 

( 

00330 
00340  _ 

A=A*EAC 

AP=-(OROOTX*RAT)»A*RTPI*ROOT4X*AP 

( 

00350 
00360 .... 

fac=eac*fac 

8 = 8*F  AC 

00370 

00380 

BP=(DROOTX-RAT)*B*RTPI2*ROOT4X*BP 

RETURN 

( 

00390 

END 

< 

• 

< 

( 

12 

11 

tO 

9 

< 

8 

< 


7 

6 
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- 00010- 
00020 
00030 
00040 

00050- 

00070 
00080.. 
00090 
00100.. 
00110 
00120- 
00125 
00130 
00140 
-00150- 
00160 
__0 01.70_ 
00180 
— 0019JL 
00200 

00210 

00220 
—002J0- 
00240 
-00250- 
00255 
_ «•£«>_ 
00270 
. 00280. 
00290 
00300 
00310 
00320_ 
00330 
00340 
003^0 
00370 
00380 
00390 
00400 
00410 
00420 
00430 
00440 
00450 


00470 

00480 

00402 

00404 

00406 

00406 


00490 

00500 

005i0 

005g0 

00530' 

005  40_ 
00550 


SUaRaUUNE-.P-R0PtAU»8E.XaJ»U<liJ.»UP^UBj4JCN£lUXUCfUll. 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-H) 

C0MM0N_/lr.LAGS4JirLAC,>1£LA6i/im£E 

COMMON  /SQNORMX  XNORM 

DM=Al*-0Q0 

OX=-AL*(X+BE) 

OXl«rAU»lXJ.+BE-> 

CALL  Oa!RY(DX,AJ.AIP.BI,BIP» 

P.I_?.3 .141592653509703238-400 

AJ*PI*( (BlP«U)-CTT*BI*UP) I v .. 

BA=P-l»UlXtAl*UP-lpAA.lP.gUll 

TT=1.0D0 

i F_(  MFL  AG  j NE  ,_i)_G0.10-5 

B1=<AJ**2)*(<(DX/AL>*(AI*»2>)-((A1Pi»«2>/AL>) 

B 2=JJB  J * * 2 )JH.(.CQX/AU  * ( B l ••2XLr.LtB  1 P?*217.AUJ 

B3=(2i0O0»Aj*BJ)*( ( (OX/AL > • ( A J *8 1 ) ) - ( (AJP*8IP)/aL) I 

BEG?-8l»B2.-B3 

CONTINUE 

I F CHF.LAG  1 . PQ , 1 > CO  TO  -4 

!E((OX.LE.OM),AND.(OX1.LE.OM>)  GO  TO  3 
C QMliNllE 

xnew=xnew*xincr 

LELLXNEM*JG1.«-X  J I GO  TQ-2 

OX=-AL*(XNEW*BE) 

LF_C OX  1LE1  DH.HR .lA.EQ.ai-GQ-TO  6 

IF(XNEU’XlNCR.LE.X)GO  TO  6 

XNE.WrXNEM-X.lNCR 

IFtMFLAG.NE.l)  GO  TO  7 

0 X«  - A.L*.IXNE  W*  BE.) 

C1b(AJ**2)*( ( < OX/AL )*< A 1**2) >»< <AIP*»2>/AL>  > 

C2*  IB  J**2 1*1  L(  DX/.AL1* (B 1 4*2 ) ) I IB IP**  2 1 /.ALU 

C3=(2.0O0»Aj»BJ>*(<(OX/AL)*(Al*BI))-t(AIP*BlP)/AL)> 

CEN MCI* C 2 ♦ C 3 

ZN0RM=-CEN-BEG 

XN  O.R  M z X N O-RM*  1N0MM 

CONTINUE 

.C  AL  L-L  ARGE  (_A.Li8ii.XNEMi-X  J .MF.jMJj.W  P.iUP  ID. 

RETURN 

CQMTJLNUE 

CALL  DAIRY(0X.AI.AIP#BI,BIP) 

I A = 1 A * l . 

WF=AJ*a:*BJ*BI 

W P.r  A J *AI  P * B J *8.1 P - 

IF<MFLaG1.E0.1)  WR l TE (NTYPE »11)  WF 

F 0 R M Alt  2P20.8) 

GO  TO  i 

XNEW=XMWtiXIMCR 

0X=-AL*<Xj*BE> 

CONTINUE 

IF(XNEm*XINCR.GT.XJ)  GO  TO  150 

XNE  W _£_XNEK*X  I NCR 

GO  TO  120 

xnew=xnew-xlncb 

CALL  PAlRY<OX,AJ,AlP,BI,BlP) 

UJ=AJ?A  L*BJ*aL 

UPJ«AJ*AlP*BJ*BJP 

I Li  MF LAJS-.  NE  .1J-  R EIURN 

C1=<AJ**2)*< ( (OX/AL) *<AI**2) >-((AJP**2)/AL>) 

C 2=JB J*»  2 ) JIUJ-DX^L L* <J.l *JL2JJ -J  (BJ P»«2)/ A L > ) 

C3=(2i0D0*AJ*BJ>*( ( (0X/AL)*(A|*BI))-<<AIP*BIP)/AL)I 


£-4. 


_00.5  60-<s 

00570  100 

00580 

00590  “ 


.HEN?Cl.*C2tC3 

FORMAT  < 7014,3) 

_ZN0RMM-.CENr.8EG 

XN0RM*XN0RM*HN0RM 

-REIURN — 

END 


00610 


rjt 


tolS  PAGE  IS  BEST  QUALITY  PRACTICABLE 


7-1- 


00010 

00020 

00030 


0004? 

00050 

00060 

00070 


00090 
00100 
00110 
00l20_ 
00130 
001.4.0_ 
00150 
00l60_ 
00170 
00180- 
00190 
.00200— 
00210 
-00220_ 
00230 
.00240  . 
00250 
00260.. 
00270 

0 0 2fl  0 
00290 
00300 
00310 
00320 
00330 
00340 
00350 
00360 
00370 
00390 
00400 
.00*10. 
00420 
00430 
00440 
00450 
00460 
00470 
00480 

00490 

00500 

00510 

00520 

00530 

00540 

00550 

00560 

00570 

00580 

00590 


SUBROUTINE  PRUL.BE,X»XJ.U»UJiUP.URJ»XNEW,XINCR^TTi 

1 X1SQ.X2S0) 

iNPLlCl.T.  OOUBLE-PREC  ISlON_CAr.H»Q-El 

COMMON  / I FLAGS/  MFLAG.MFLAG1.NTYPE 

COMMON  /JC/-JCOUNI 

X1SQ»0.0O0 

X2SQ=0.0O0 

DM=11.0D0 

__IA  = 0 

OX=-AL*(X*OE) 

DX1»-AL*(XJ*BEJ 

CALL  OaIRY(OX.AI.AIP.BI»BIP) 

PI =3^141592653589783238400 

AJ=PI«< (BIP*U)-(TT»BI*UP>) 

B J=RL«  L(XL»  AJJLUP1  ■ LAIEgU-U 

TT=1.0O0 

I F1MFL AG . NE_D_G0-I0-5 

B1=(AJ*»2)»< (<0X/AL)*<AI**2))-< (AIP*«2)/AL>) 

B2=(BJ»»2)«4.(.C0X^AU»CBI«*2))-C(BIP.*<i21/ALl»- 

B3=(2,0O0»AJ«BJ)*<< (0X/AL)»{AJ*BI))-((AIP*BIP)/AL)) 

BEG=-Bl-B2-a3 

CONTINUE 

1EIMFLAG1  iEG^LI GO_Ifl— 4 

IF((0X.LE.0M).AN0,<0X1.LE,0M))  GO  TO  3 

CONTINUE 

XNEWeXNEW-XINCR 

I F.(  X N E.  w , LI  r*  JT  _ G C._HL_2 

OX=-AL*(XNEW*BE) 

I F (DX , LE . DM1_GQ--IQ— 6 

XNEW=XNEW*XINCR 

I F < MFL  AG . NE,  D-  GO  ..IO_Z 

OX=-AL*(XNEW*BE) 

Cl=  I AJ«.*2 ) * UXOX/AL  > * ( Al**2).Lr.UA  IP**2  i/ALU 

C2= < B J*»2 ) a ( { <0X/AL)MBI**2) )"{ <BIP**2)/AL> I 

C3  = (2i  0Q0»A  J*BJJ«14-(JDX/AL-).«tAL*B.l)lr.(  lAlPtfllP.I^ALU 

CEN=C1+C2*C3 

XlS.O  = CEN.t.BEG 

CONTINUE 

L FORMA  It.'— PR.'.  a5015^1  J __ . , 

CALL  iGtAL.BE.XNEW.XJ.WF.UJ.WP.UPJ.XNEW.XINCR.TT; 

1 X2SQ<  T2SQJ 

RETURN 

CONTINUE 

CALL  0AIRY(0X,AI.AIP.BI»8IP) 

WF_=AJ*A  J*BJ*BI 

WP=AJ*AlPaBj*BlP 

J F ( M F_  LAG  li  E Q «J,  J _WRl_lE.<iili.PE_i  Hi— ME 

FORMAT ( 1020 i 8 ) 

GO  _TQ__i 

XNEW=XNEW*XINCR 

DX  = -A  L • ( X J ♦ BEJ 

CALL  OaIRYIOX.AJ .AlP.Bl ,BIP) 

UJ=AJ»AI*BJ*BI 

UPJ=AJ*AIP*BJ*BIP 

iflmflag.ne.ilreturn 

C1»(AJ**2).( ( (0X/AL)*<AI**2))-<(AIP**2)/AL>) 

C2s<BJ**2j_«U  ( 0 X / A L ) •_<  B 1*  * 2 1 l-J  < 0 1 Pi « 2J/-A L>I 

C3=(2,0D0*AJ*BJ)«(((OX/AL)»(AI*BI>)»((AIP*BIP)/AL)I 

CEN»Cl*C2aC3 

X1S0=CEN*BEG 


00810- 

00620 


. RETURfi 
ENO 


f 
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..  00010. 

.SUBROUTINE  START  CAL.  BE , X.  XJ»  Ua  UP.iXNEW«-)t!NCR.-SW..tTL 

00020 

IMPLICIT  DOUBLE  PRECISION  (A-H.0-2) 

00030 

COMMON  /IFLAGS/  .MFLAG . MFLAG1. NT YPE  . 

00040 

COMMON  /SQNORM/  XNORM 

00050.. 

nXr-AL*<X*BF ) 

00060 

dxi=-al**<xj^be> 

00070 

irnOX.GT.15. 000).  OR.  (0X1.  Cl.  15..0O0D  GO  TO  4 

00080 

CALL  OaIRY<OX.AJ,A!P.BI*BIP) 

00090 

a=sw*Bi*tt, 

00100 

B=-At*SW»TT 

00110  .. 

A J= A ... 

00120 

BJ=8 

...  00130  .. 

IF (MFLAG »NE ..1)  GO  TO  5 . 

00140 

B1  = <AJ**2)*(  ( (DX/AD»(Al**2))-(  {AIP**2)/AL>) 

...  00150 

..  B 2 = ( B J • »2  > * LL  ( Q X / ALI  • IB  I* « 2 ( 8.LB*  *2 >./  AL ) 1 

00160 

B3=(2,0D0»AJ*BJ)»(<<OX/AL)*<AI*BI > )-( (AIP*BIP)/AL)> 

... .00170  . 

BEG=-@L-0  2.-33. 

00180 

5 

CONTINUE 

. . 00190  . 

OXs.-AL*tXj.*BEl 

00200 

CALL  OaIRY(OX»AI.AIP.BI,BJP) 

<?0?1_0. 

_ u=a»ai*b*bi 

00220 

Upz a»a i p*b*b I p 

00230 

IF  (MFLAG iNEUL -RETURN 

00240 

2A*AJ**2 

.00250  . 

ZB= (OX/ALI 

00260 

2CaAI«»2 

00270 

20=<AlP**2)/AL 

00280 

C1=ZA*( (2B»2C)~ZD) 

00290 

Cl* 1 A J* •2 1* 11  (DX/AL)  «.£.A.I*.*2.)  )_-!( AIP*»2 ) /AL ) ) 

00300 

C2=(BJ**2><*<  ( (0X/AL)*(BI»*2) )-(  181P»«2>/AL> ) 

. 00310 

. ...  C3=  ( 2 « 0Q0»Aj«BJ ) » LU 0X/AL1  * UUB LI  >-t<ALP*B IPJ / ALU 

00320 

CEN=C1*C2*C3 

20330. 

2NCRM  = -CE.\'-BEG 

00340 

XN0RM=XN0RM*2N0RM 

00350 

1 

XNEW=XKEW*X1NCB  

00360 

IF(XNEw.GT.XJ)  GO  to  3 

00370 

DX=-AL*<XNEW*BE) 

00380 

CALL  DaIRY(DX,AI.A1P.BI,BIP) 

00390 

WF  = A » A I *B*B I 

00400 

IF((MFlAG1.EQ,1).AND. (XNEW.GT, 0,000))  WRITE(NTYPE.ll) 

00410 

1 WF 

00420 

Up=A«A I P*B*B I P 

00430 

u 

FORMAT(2Q20..81_ 

00440 

GO  TO  1 

00450 

3 

XNEW=XNEW-XINCR 

00460 

return 

00470 

4 

CONTINUE 

00400 

c 

OXs-ALo(X+BE) 

00490 

c ... 

CALL  EXTReM(DX.A2.AZP.BZiB2P»H>  .. 

00500 

XNEWsX 

00510 

T=1.0O0 

00520 

c 

A*BZ*Sw 

00530 

c 

B=-AZ*SW 

00540 

c 

U1=A*AZ*B»B2 

00550 

c 

UP1=A»AZP*B*B2P 

00555 

U1*0 

00556 

UPl=-TT/3, 141592653589783200 

00555 

c 

TYPE  U.X.OX.UliUPl 

00500 

CALL  LARGE(AL»BE.X.XJ*U1.U»UP1,UPiXNEW,XINCR.T) 

00590 

RETURN 

ENO — 


00600.. 


THIS  PACE  IS  BEST  QUALITY  PRACTICABLE 
FROM  COPY  FUtraUSHED  TO  DDC  


00010 

00020 

0003® 

— SUBROUTINE  LARGE(AL#BE.X»XJ*U«UJ*UP# UP J»  XNEM*  XlNCR«  T-T 1 

IMPLICIT  DOUBLE  PRECISION  U-H.0-Z) 

COMMON  ZlFLAGS/  Mf LAGiMFL AG1 1 NTTPE 

00040 

00050 

common  /sqnorm/  xnorm 
□ Ms  11 < 000 

00060 

00070 

U INI T"U 

...  . SUMQ  = 0.0O0_  ...  _ ....  _ _ 

00000 

00090 

SUMEs0.0O0 

IAs0  

00100 
00110 

PI =3, 1415926535889783238400 

OX  = -AL«.(X*BD 

00120 

00130 

1 . 

CALL  ExTREM(0X,AZ,AZP.BZ.B2P,2> 

..  XNFW=XNEW*X?,NCR  . ...  _ . ..  .. 

00135 

00140  _ 

IF(XNEw.LT.X)  GO  TO  1 

IF(XNEw.GT.XJ)  GO  TO  2 

0015® 
00160  .... 

DX=-AL*<XNEW*8E) 

IF  (DXtGT  i.OM)  GO  TO  4 _ 

00170 

00100 

xnew=Xnew-xincr 

_ . IF (MFLAGtNE.l)  GO  TO.  5 ... 

00190 

00200 

1F( IDF ,EO,0)  SUME=SUME-(4,0D0»U1»U1) 

IF( IDF.NE.0)  SUM 0= SUMO- (2.0D0*U1*U1) 

00210 
00?20  ... 

SENO= <UINIT**2)*(U1**2) 

2NORMs  t X I NCR/3 . 0D0  ) #_t  SUME* SUMO* SEND ) 

00230 

00240 

5 

XN0RMSXNGRM*2N0RM 

CONTINUE 

00250 

00260 

CALL  PR0P< AL.BE,XNEW.XJ,U»UJ,UP|UPJ,XNEW»X1NCR.TT) 

RETURN 

00270 

0028® 

4 

CONTINUE 

CALL  ExTREM(DX»A22.A2P2iB22»B2P2»H2> 

- 

01*290 

00300 

DEX=DEx?<2-22) 

C11  = PI*C  <A22«BHP*0EX>-<  <922*A2P)/0EX> ) 

00310 

00320 

C12=<PI*TT)*( ( (B22»A2)/0EX)-(A22*8Hm0EX) J 
C2ls(Pl/TT)*((A2P2*B2P»DEX>-UB2P2*A2P)/DEXn 

003  30 
00340 

C22=PI«(<{82P2»A2)/0EX)-(A2P2*B2*0EXW 

U1=C11*U*C12#UP 

00350 

00360 

11 

IF(MFLAGl.EO.l)  WRITE(NTYPE.ll)  U1 

FORMAT ( 2020 . 8 ) 

00370 

00380 

IF(MFLAG.NE.l)  GO  TO  3 

I A = 1 A*1 

00390 

00400 

I H* I A/2 

. . iDF=«  1A-2*III>  . 

00410 

00420 

!F( 1DF.EO.0)  SUME=3UME*C4,0O0*U1»U1) 

IF(IOF.NE.0)  SUMOrSUMO* (2 • 0D0*U1*U1 1 

. . 

00430 

00440 

3 

CONTINUE 

UlP=C2l«U*C22»UP 

00450 

00460 

UlP3UlP*TT 

TT»1,0O0 

00470 

00400 

OXs-ALa(XnEW«BE) 

A2=AZ2 

00490 

00500 

AZP* AZP2 

BZsB22 

12 

11 

00510 

00520 

8ZP=BZP2 

2 = Z2 

U 

00530 

00540 

U*Ul 

UPsUlP 

e 

00550 

00560 

100 

FORMA  T(4Dl5,6) 

GO  TO  1 

a 

9 

00570 

00580 

2 

continue 

xnew=xnew-xincr 

00590 

DXs-AL*<Xj*BE) 

1 


3 


10- 1 


00010.-  

00020 

00030 

00040 

.00050 

00051 

00060 

000Q0 

00090— 

00110 

00130 

00140 

00150 

00160 

00170 

00180 

00190. 1 ... 

00200 

00210 

00220 

00230  - 

00240 

00250. 

00260 

00265  

00270 

00280  14 

00285 

00290 

003n)0  t5 
00305  20 

003i0 
00320 

00330  10C 

00340  _ 5 

00350 

00360 

00370 

00380  4 

00390 

00400 

00410 

00420 

00430 

00440 

00450 

00460  11 

00470 
00400 
00500 

00510  502 

00520 

00530 

00610 

00620 

00628  " 30 
00630  __  _ 
00640 

00650  

00660  4*. 
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SUBROUT  INE-LGCAL*8E»X,.XJ.U.UJ. UP.  UP  J J(NE«,XlNCR,-i 
X1S0.X2SQ) 

IMPLICIT  double  PRECISION  (A-H.O-H) 

COMMON  / I FLAGS/  MFL AG , MFL AG1 . NT YPE 

COMMON  ./ JCZ-.JCOUNT 

DIMENSION  W ( 250 ) » WP  ( 250 ) 

...  OM=11.0D0 

X1SO=0.0O0 

X2SQ=0«0D0 

I A = 0 

UINIT'U 

SUMO=0.0D0 

SUME=0.0Q0 

PI =3. 1415926535889783238400 

._  0X=-AL«(X*8E) 

CALL  ExTREM(DX»AZ.AZPiB2»BZP»Z) 

XNEW5XNE.W.-rXLHCR 

FCT=1.0O0 

IF.tXNEW.LT.XJJ  GO_I0_2 

OX=-AL*(XNEW*BE) 

_ iFlOX.GT.DM.OR.  IA,LT..1)_G0  .tU-4 

XNEW=XNEW*XINCR 

„ IF  1MFLAG  «NE,_L)  _GQ_IQ.-5 

1M=IA-1 

...  IF(IM,LT,1)  GO  TO  .20 

00  14  K=liIM,2 

IF(OABS(W(K) ) .GT.l.D-18).  SUM05SUMO+ (2 . 0O0»W  IK ) «*W  < 
IFdM.LT.  2)  GO  TO  20 

_ DO  15  K?2,  IMj2 

IF  (OABS  ( W (K) ) .GT.l.D-18)  SUME  = SUME*  ( 4 . 0D3*W  { ) ow  ( 

CONTINUE  ....  

SEN0=(UINIT*»2)*W( IA>«»2 

. . XlSO-<XlNCR/3,0O0.)«tSUME*SUMQ*SLNO) 

FORMA  T ( 5014 , 3 ) 

CONTINUE , 

CALL  PR  I Al.BE,XNEW.XJ*UiUJ.UP»UPJ*XNEW.XINCR,TT, 

. X2SQ.Y2SQ) 

return 

CONTINUE _ 

CALL  EXTREM<0X.AZ2.AZP2,BZ2,BZP2»Z2> 

_ OEX  = OEXP  ( JL-Z2 ) 

Cll  = PI»dAZ2*BZP«DEX)-(  <BZ2oAHP)/DEX)  ) 

_ .C12=<P.1*TJU<  ( <BZ26AZl/DEX)r(AZ26BH»DEX)l 

C21=(Pl/TT)*( (AZP2»BHP«0EX)-( (BZP2oAZP)/0EX> ) 

....  C 2 2 = P ! * <J  ( B Z p 2 * A Z } / 0 E X ).r 1 A 2 P 2*B  1*0  EX ) 1 

U1=C11«U*C12»UP 

...  FORMAT  (1023. 8) 

U1PMC21*U*C22*UP)*TT 

TI  = 1 ,000 

1 A = I A ♦ 1 

FORMAT!  (3j 

W(IA)=U1 

WP ( I A ) ?UiP 

U=W(1 A) 

UPrWPUA) 

F0RMAT( 13,01 

!F<MFLAGl,EQ.l > MR  t TE ( NTYPE .11)  W( 1A  I 

OX=-AL»<XNEW*BE) 

AH* AZ2 

AZP-AZP2 


( 


\ 
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TOO*!  COPY  FURNISHED  10  UDQ  

?0670 

02  = 822  - _ . _ ..  . . ....  . 

00680 

00690 

BZP»BZP2 

2 = 22 

\ 

00700 

00710 

2 

CO  TO  1 

CONTINUE  . - . ...  

00720 

00725 

IF(MFUAG.NE.l)  GO  TO  3 

I F C IAtLT.l)  CO  TO  3 - - 

{ 

00730 

00740 

i6 

00  16  KB1 » I A » 2 

1F(DABS(W(K) ) ,GT. 1.0-18)  ,SUM0  = SUM0*2 . 0D0*W (K ) ••2.r  . 

j 

< 

00745 

00750 

I F ( IA.LT.2)  CO  TO  3 

DO  17  K = ?,  1 a i 2 ...  - . . _ . . 

00760 

00770 

17 

3 

IF(0ABS<W(K)  ) .CT.l.D-18)  SUME  5SIJME  + 4 . 0O0«W  < K > **2 

CONTINUE  _ ._  . . ... . ... 

00780 

30790 

xnew=xnew+xincr 

OX=*AL*(Xj*BE) 

/ 

00795 

00798 

26 

I F ( OX  t L T i 5 > TYPE  28, OX 

FORMAT  ( i WARNING*.  IN  LG  ARG  OF  EXTREM  IS  . ' , D I 

00600 

30610 

CALL  ExTREM(0X.A22,AZP2,B22,8HP2,22) 

DEX=DEXP(Z-Z2)  - . . 

1 

00820 

30330 

CllaPI»(  I AZ2°B2P'*DEX  ) - ( ( 8Z2* A 2P ) /OEX  ) > 

C123(Pl«TT)»({(BZ2«AZ>/0EX>-(AZ2»BZ*0EX)) ...  . 

{ 

00340 
00  850 

C21=<PI/TT)»((AZP2«BHPoDEX)'((BHP2»AHP)/0EX)) 

C22  = Pl»(  ( (B?P2*AZ)/QEX)'-CAZP2«B2»0EX1).  . ...  

00860 

00870 

UJ=C11»U*C12»UP 

UPJ=(C21<*U*C22»UP)*TT  .......... 

{ 

00380 

00930 

IF(MFLAG.NE.l)  RETURN 

SEN0=(UINIT*»2)*IUJ*°2)  ...  . ......................... 

00940 

00950 

X1SQ=<XINCR/3,0O0)*(SUME*SUMO*SEND) 

RETURN  ....  . 

00  6 0 

END 

' 

( 
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< 

12 

11 

10 

1 

